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INTRODUCTION 


Recently a report [1] was prepared describing the theoretical under- 
pinnings of a vortex wake method of analysis for prediction of the aerodynamic 
performance of horizontal axis wind turbines. By use of this method, rotors 
having any number of blades, which may be arbitarily shaped and twisted, can 
be analysed. Furthermore, it has been found that a computer code entitled 
VORTEX which implements the vortex wake method of analysis can obtain answers 
to problems of interest in computing times that are comparable to those 
obtained from more elementary methods e.g. the blade element -momentum theory 
used in the well known PROP computer code [2] or in the WIND II code [3], 

In the vortex wake method, the induced velocity is directly obtained by 
integration of the Biot-Savart law. To implement this integration, it was 
assumed that a discrete number of vortex filaments trail from the rotor blade. 
These filaments extend infinitely far downstream and have a constant diameter 
helical shape. It was also assumed that the entire helical vortex system 
produced by the rotating blades travels downstream with a constant velocity 
equal to the value of the rotor disk. Lifting line theory was used to 
represent the blades and aerodynamic effects were incorporated by use of 
empirical airfoil lift and drag curves. 

Aerodynamic performance predictions using the vortex wake method of 
analysis were compared to experimental data for four different rotors in [4]. 
In general, favorable agreement between measured and predicted rotor power was 
obtained. The method has also been recently extended to include performance 
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prediction of tip-controlled wind turbines [b], [6] and [7J. Moreover, a 
vortex wake model was found to provide better simulation of an experimentally 
observed post-peak power plateau than did a blade element -momentum model [8], 

A significant drawback of the calculation procedure needed to implement 
the vortex wake method of analysis is that it is unavoidably quite involved. 

The reasons for this complexity are better understood once the overall pro- 
cedure has been outlined. Essentials of the procedure were given in [1], 
however, certain important numerical details were omitted in order to avoid 
complicating the presentation of the theory. In addition, no discussion of the 
computer code that evolved from that study was presented. As a consequence, 
the current report is written with two principal objectives in mind: 

(1) to supply missing details in the computational procedure, and 

(2) to describe the computer program VORTEX developed to implement the 
vortex wake theory. 

In order to accomplish these objectives, both the numerical procedure and 
the computer program will be fully discussed in the following. Details of the 
integration and interpolation schemes that were used will be presented along 
with a method for treating numerical singularities that occur in some of the 
governing integral expressions. All program variables and subroutines will be 
defined and input/output parameters will be described. Finally, program 
implementation will be illustrated by the presentation of an example problem. 
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THE NUMERICAL PROCEDURE 


For calculation purposes each blade of the wind turbine is thought to be 
divided into M-1 spanwise sections. These subdivisions establish M calcula- 
tion positions along each blade. Although the subdividing is arbitrary, it 
must be understood that a large number of subdivisions can lead to excessive 
computational times without necessarily improving overall accuracy. In most 
problems solved to date, 9 subdivisions were found to provide good balance 
between accuracy and computer run times. It should also be mentioned that the 
subdivisions need not be of equal size e.g. , more sections may be located near 
the blade tip or in regions of high gradients. To simplify the computations, 
it is presumed that each of the N rotor blades of the wind turbine has been 
subdivided into the same number and distribution of parts. 

In the vortex wake method of analysis, it is assumed that each rotor 
blade can be represented by a single bound vortex i.e., the lifting line. And 
because no vortex line may end abruptly within the fluid, it is assumed that 

there is a system of vortices that trails each rotor blade. Accordingly, it 

is appropriate to think of a vortex filament emanating from end points of a 
blade subdivision. The collection of these trailing vortex filaments form a 
system of vortex sheets having a helical geometry due to the blade rotation. 
These vortex sheets extend infinitely far downstream of the rotor and are 
assumed to have a constant diameter and constant pitch; this is known as the 

rigid wake assumption. In point of fact, the wake is not rigid but expands 

radially in the downstream direction. 
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In performing the calculations, the influence that each trailing vortex 
filament has on the flowfield in the vicinity of each bound vortex must be 
determined. The combined influence is called the induced velocity or the down- 
wash. The downwash can be found by direct integration of , the Biot-Savart 
law [9]. This integration is made somewhat difficult by the fact that it must 
be performed over the entire length of every helical vortex filament. The 
downwash is subsequently used to calculate the induced angle of attack distri- 
bution which in turn may be used to evaluate various performance factors, e.g., 
rotor power, rotor thrust, etc. To accomplish all of this, the circulation 
distribution along the blade, r(c), where 5 is the nondimensional spanwise 
position dimension, must be known. As will be shown in the following, the 
particular form of the distribution evolves from the calculation procedure 
having initially been assumed in the form of a truncated Fourier sine series: 


ST n , ^ ~ 5hub— I 

P’Z f-ThuTJ 

m=l 


( 1 ) 


It should be noticed that this distribution has been constructed so that the 
circulation at both the blade tip (5=1) and the blade hub (5=?hub) vanishes. 

The lift and drag coefficients for each blade section must also be known. 
Generally, this information is supplied in the form of curve-fitted wind tunnel 
airfoil data. 

In order to have an understanding of the calculation procedure, the 
following step-by-step outline and description is presented. 
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Step 1 Select M stations along the blade span. 


Each calculation site on the blade is designated by a 5 ' value. As men- 
tioned, generally nine (9) stations have been used in the computations. For the 
numerical integration, a separate nodal system of blade positions is required. 
These locations will be designated as Obviously, there are many more 5 loca- 
tions than 5 ' locations. 

Step 2 Determine the circulation distribution along each blade (first pass). 

In the calculation procedure, the circulation distribution was determined in 
one of two ways depending on whether it was the first pass through the calculation 
procedure. In the first pass, it was found helpful to write the effective angle 
of attack for each blade section, {oie)^,, as a function of sectional lift coeffi- 
cient, (Cl)[^. This expression when combined with: (a) the Kutta-Joukowski 

Theorem [9], (b) an expression for the induced angle of attack and (c) the circu- 
lation distribution produces a rather complicated expression {see equation (2-67) 
in [1]} which can be written in compact form as 
M 

^ fU')A„ = g(c') (2) 

m=l 

In this equation, both f and g are functions that involve many parameters and 
the are the coefficients of the circulation distribution in equation (1). By 
applying equation (2) at each 5 ' location, a set of M equations in M unknowns is 
obtained. Solution of that set gives the A^, which permits the circulation to be 
determi ned. 

Step 3 Calculation of the induced angle of attack. 

At each station, 5 ', the induced angle of attack, computed 

from the circulation distribution. 
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Step 4 Calculation of the effective angle of attack. 


From the relation between geometric, effective and induced angles of 
attack, new values of can be evaluated. 

Step 5 Calculations of the lift coefficient. 

Using two-dimensional airfoil data and the values of the effective angle of 
attack, new values of (C|_)fp can be determined. If these agree with the previous 
(CL)m» the iteration procedure is terminated. If thoy do not agree, the process 
continues to the next step. 


Step 6 Determination of the circulation distribution. 

The Kutta-Joukowski theorem applied to each blade section may be written 

^m “ 2 

Values of (CL)fn from step 5 permit the to be determined. In turn, the A^, in 
equation (1) may be found from the following matrix equation 


[C] {A„,} = {r} 


( 4 ) 


where 


Cij “ sin 


, ^ / 1 ■ ?hub 

(J^) (-r )| 

^ ■ ?hub 


Once this calculation is completed, the program is directed to return to step 3. 


6 



NUMERICAL ANALYSIS 


Because the computer program must perform interpolations and integrations, 
solve systems of equations and resolve numerical singularities that arise in 
certain integrals, it is appropriate that this numerical work be described. 

Interpolation 

In the program, a cubic spline interpolation technique [10] was used in 
which the second derivative at each end of the interpolated data set is assumed 
to be a linear extrapolation of the value at the two adjacent points. This 
interpolation was found to be necessary to calculate the induction factor in the 
neighborhood of the singularity (see singularity treatment below). 


If a cubic polynomial is defined as 

f(x) = + a^x + a^x2 + a^x3 

then between points X|^ and X|^+i, the second derivative has the value 


, .. /Xk+1-X\ .. /X-Xk 

(X) = f k ) + f k+1 

\ 5k / \ 5k 


( 5 ) 


where 6k = Xk+i - Xk. Integrating equation (5) twice yields 


II 

~(XkH - X)’ 

M 

(X - Xk)^ 

f k 

_ 6 6k _ 

+ f k+1 

__ 6 6k _ 


+ C X + C (6) 
1 2 
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The integration constants and are obtained from the fact that f(X) passes 
through X|^ and X|^+i. Thus, 


(^k+1 - ^k) (f"k+l - f"k)^k 



(^k^k+1 “ ^k+l^k) k^k+1 * ^ k+l^k)^k 

6k 6 


Inserting these into equation (6) produces 


f(X) 


f"k(Xk+l - 

6 6k 


f"ktl(x - X|()3 

6 6k 


+ (Xk+1 - X) -- 
^k 


f"k^k 


+ (X 



^"k+l^k \ 

....... j 


( 7 ) 


In this equation, only the second derivatives f "k and f k+1 unknown. 
However, the value of the second derivative at all points can be obtained by 
equating the slope of two neighboring subdomains at both ends of the inter- 
polated data set. This results in m-2 equations for m points. Therefore, two 
additional equations are required. They may be obtained by linearly 
extrapolating the value of the second derivative at the two points adjacent to 
both ends of the interpolated domain i.e. , at points 1 and m. This results in 
the following expression for each internal point: 
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f"k-l «k-l / \-l ^ 

6 -- * ' \ T 

and the following pair at the end points 



^ k+1 ^k _ ^k+1 ■ ^k 

6 6k 


^k - ^k-1 
^k-1 


( 8 ) 


(9) 


( 10 ) 


This system of linear equations can be solved for f"k» k = 1 m. Thus, 

the interpolating function is completely determined as 

f(X) = Ai_k(X|,ti - X)3 + A2,|((X - X|,)3 t A3,k(Xk+i - X) 

t A4 _h(X - X|() (11) 


Where the are the coefficients of the x values as displayed in equation (7). 
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Integration 


The definite integrals that occur in the analysis are handled by making 
use of the spline interpolation formula, equation (11) above. The integration 
formula that is produced is 


/' 


Ai 1/ , ^2 k 

f(x)dx = (Xk+i - X)** + (X 


Xk)** + 


Oik+i - X)" + (X - 


( 12 ) 


The semi -i nf 1 nite integral that occurs in the equation for the induction 
factor was determined by using a 24 point Gauss-Laguerre formula, i.e.. 



f(X)e“XdX 


0 



(13) 


The weighing factors, Wj, can be found in the literature, [11]. 

Solution of the governing system of equations: 

The linear system of equations, equation (4), for the circulation coeffic- 
ients IS solved simultaneously by employing the Gauss elimination method. 
Because this method is rather elementary, it is not presented here. Those 
unfamiliar with the technique should consult a standard numerical analysis 
reference, e.g., [12]. 
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Treatment of the numerical singularities. 


In [1], it IS shown that the total induced velocity at any point ^ on the 
blade lifting line due to all trailing helical vortices from all blades is 



/ - dr 

^tip --dc N /•«> 

Y / Ni?' + 

Z- / d‘3/2"d 

■^nub ><=1^0 1 


^2^0 

i/z 


do 


(14) 


where Nj 


N2 


Dl 


D2 


5[? - 5'cos(e + 0|<)] 

[-?' + K cos(0 + 0k) + 0 ? sin(0 + 0|<)] ^ 

K 

+ 5 '^ - 2^5'cOS(0 + 0 k)+ --9^ 


h 


^0 ■ Wn 


Q + 


COS(f) 


s i n<|) ' 


r 


I 


Careful examination of equation (14) reveals that a numerical problem can 
develop. In particular, if the parameter Dj becomes zero, as it can when 
0 = ©k = 0 and C = 5 ' occur simultaneously, the value of the induced velocity at 
that point on the blade becomes infinite. The presence of these points of 
singularity prevents a direct solution of equation (14) in the neighborhood of 
the poi nts. 
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To overcome this numerical difficulty, a theory developed by Monya [13] 
was utilized. This theory is based on the fact that if the differential normal 
induced velocity [the differential form of equation (14)] was written for a 
straight trailing vortex (as opposed to the actual helical vortex), this quan- 
tity would become unbounded in exactly the same manner as the quantity for a 
helical vortex. Physically this results from the fact that the major contribu- 
ting factor which causes equation (14) to become unbounded is due to the segment 
of the trailing vortex that is very close to K' i.e., very close to the trailing 
edge at c'. Clearly, it is assumed that curvature effects of the helical vortex 
in the vicinity of the singularities are not important. Figure 1 is a graphical 
representation of this concept. Since the differential normal induced velocity 
for both the helical vortex [say (dWp)H] and the straight vortex [say (dWn)s] 
approach infinity at the same rate , then their ratio, which Moriya called the 
induction factor I, i.e.. 


(dWn)s 

tends to unity at each singularity. 


Combining equations (14) and (15) along with an expression for (dWn)s, [9]. 
yields 




5tip dr 


d? I 

4ttR U - 


(16) 


^hub 
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where 


I = 1(5, 5 ', Xq) 



Nl?' + N2X0 

d3‘/T"i72" 
1 2 


( 17 ) 


Experience with the numerical integration of the induction factor using an 
integration interval spacing equal to ten times the number of blade segments 
minus one, times the blade radius revealed a problem. It has been found that for 
small values of the wind-tip speed ration, Xq, the numerical value of I tends to 
oscillate both slightly before and slightly after the location of the singular 
point. This behavior is diagramed in Fig. 2a. To remove this oscillation in a 
sensible manner, it was decided to integrate equation (17) to within a small 
distance of the singularity (denoted by I = 1) and then to use the last four I 
values to fit a cubic spline through the point I = 1 as diagrammed in Fig. 2b. 

The distance from the singularity needs to be large enough so not to contain any 
of the unwanted oscillation yet not so large as to result in a loss of accuracy. 
To be sure, the selection of the distances on either side of the singularity 
requires some experimentation. 

It should be noted that Moriya did not report encountering any oscillation 
of the type found in this study. This is not surprising since his work concen- 
trated at much larger values of Xg where no such oscillation has been observed. 
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FIG. 2a OSCILLATION IN THE INDUCTION FACTOR INTEGRATION 



FIG. 2b REMOVAL OF OSCILLATION IN THE 
INDUCTION FACTOR INTEGRATION 
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THE COMPUTER PROGRAM 


A computer program to implement the numerical procedure described in the 
previous section has been written in FORTRAN IV. Figure 3 displays a flowchart 
of that program. It can be seen that the program is divided into five primary 
subprograms. These subprograms in turn call 11 other secondary subprograms. 

In the following each of the primary subprograms will be described in order 
of their appearance in the program. This description will consist of: (a) a 
brief outline statement of the function of the subprogram, (b) a listing of the 
subprogram arguments and local variables, (c) a list of all secondary subprograms 
called and finally (d) a collection of any important notes relevant to the 
particular subprogram. 

A similar description of the secondary subprograms, listed alphabetically, 
will also be provided. 

The Primary Subprograms 
ASSGN 


Function: The purpose of this subprogram is to process input and 

output data and to calculate necessary program constants. 

Arguments: All arguments are transmitted by COMMON except 

OMEGA = rotational velocity (q), RPM 
VEL = wind speed (V), mps 
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^ cases ^ 
Computed 


FIG. 3 FLOW DIAGRAM FROM THE COMPUTER PROGRAM VORTEX 
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SI = coning angle (\|<), degree 

RDC = radians to degrees, conversion constant 

NB = number of blades 

PI = t: 

R = parameter defined by Rcosij^ 

SKPl = interval ahead of the singularity where the 
induction factor is interpolated 
SKP2 = interval after the singularity where the 
induction factor is interpolated 

Local Variables: 

BT = twist angle (e), degree 

Secondary subprograms called: 

SKIP 

Notes: To avoid recalculation of constants, all required constants have 

been calculated in this subroutine and transmitted to the main 
program. Also this subroutine calculates the parameters SKPl and 
SKP2 if ISKP is set equal to 1; otherwise, the interval for inter- 
polation around the singularity is read into the program. 
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CINTL 


Function: The purpose of this subprogram is to establish the matrix 

equation to be solved for initializing the Fourier sine series 
coefficients. The coefficient matrix is constructed based on a 
linear lift line approximation. 

Arguments: All arguments are transmitted by COMMON except 

NP = number of points used in the spline interpolation and inte- 
gration. This number must be at least 4 * N so that there will 
always be at least 4 points to obtain spline constants 
(subroutine SPLCOE will fail with less than 4 points). 

ZETA = non-dimensional radial position along the blade 

ROC = radians to degrees, a conversion constant 

ALG = geometric angle of attack, «g 

TIGRL = matrix containing Fourier series terms 

CA = coefficient matrix 

RAS = known vector 

Local variables: 

DELI = a small number, 10"® 

YI = array containing the induction factors 

Secondary Subprograms called: 

TINGRL 
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ALPHAE 


Function: The purpose of this subprogram is the calculation of the 

effective angle of attack, ag, induced angle of attack, , 
and circulation distribution, r, along the blade using an 
iteration method. 

Arguments: All arguments transmitted by COMMON except 

CA = coefficient matrix in matrix equation 
RAS = known vector in matrix equation 
TIGRL = matrix containing Fourier series terms 
R = parameter defined by Rcosij; 

ALG = geometric angle of attack 

ALE = effective angle of attack 

ALI = induced angle of attack 

GAMM = circulation distribution function 

Local variables: 

VR = undisturbed resultant velocity 
ARG = argument of sine and cosine functions 
MAXITR = maximum number of iterations 
NITER = index of iteration loop 

STIGR = parameter containing the integral part of the induced velocity. 

Secondary Subprograms called: 

CD230 and GAUSS 

Notes: This subprogram contains the following error message - "convergence 

was not obtained within the specified number of iterations". 
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CALPWR 


Function: The purpose of the subprogram is the computation of performance 

parameters, axial force, torque, power etc. 

Arguments: All arguments are transmitted by COMMON except 

ALE = effective angle of attack 
R = parameter defined by Rcos\|) 

Vel = wind velocity 

AFRCE = distribution of axial force on the blade 

TRQ = distribution of torque along the blade 

SUMF = total axial force 

SUMQ = total torque 

CF = axial force coefficient 

CP = power coefficient 

RPWR = rotor power 

APWR = alternator power 

XRA = inverse of tip speed ratio 

Local variables: 

PHI = inflow angle (angle between resultant velocity and the axis of 
wi nd turbi ne). 

VR = undisturbed velocity 

CY = parameter defined by Cl coS(i) + Cp sin<j> 

CX = parameter defined by Cl sin<(> - Cp cos (j) 

RHB = hub radius 

Secondary Subprograms called: 

CD230, CD44, SPLINT 
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OUTPUT 


Function: The purpose of this subprogram is to write out the performance 

parameters. 

Arguments: All arguments are transmitted by COMMON except 

ALE = effective angle of attack 
ALI = induced angle of attack 
OMEGA = rotational velocity (n), RPM 
R = parameter defined by Rcosip 
VEL = wind speed 

AFRCE = axial force distribution 

TRQ = torque distribution 

SUMF = total axial force 

SUMQ = total torque 

CF = axial force coefficient 

CP = power coefficient 

RPWR = rotor power 

APWR = alternate power 

XRA = inverse of tip speed ratio 

Secondary Subprograms called: 

None 
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The following describes the supporting secondary subroutines of the main program 


AUX 

Function: The purpose of this subroutine is the calculation of the 

integrand of the semi -inf mite integral. 

Arguments: All arguments are transmitted by COMMON except 

THET = independent variable of the cylindrical coordinate system 
FX = integrand of the integral 

Secondary subroutines called: 

None 


CALCI 

Function: The purpose of this subroutine is the calculation of the 

induction factor. 

Arguments: All arguments are transmitted by COMMON except 

XI = induction factor 
NOPT = control parameter 

(when NOPT = 2, calculation is performed for second blade 
NOPT = 1, calculation is performed for first rotor.) 

Secondary Subroutines called: 

GLQUD 
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CD230 


Function: The purpose of this subroutine is the calculation of the lift 

and drag coefficients using curve-fitted empirical data for NACA 
230XX airfoils. 

Arguments: All arguments are transmitted by COMMON except 

IREN = control parameter 

(when IREN = 0, no Reynolds number effect is considered, 
when IREN = 1, airfoil characteristics are corrected for the 
Reynolds number) 

ALPHA = effective angle of attack 
CL = lift coefficient 
CD = drag coefficient 
X = non-dimensional radius 

W = relative velocity 
Local variables: 

All airfoil parameters are referred to a smooth airfoil of 18% thick- 
ness to chord ratio at a Reynolds number of 3 x 10®. Some of the 
parameters are shown schematically in Fig. 4. 
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FIG. 4a LIFT 0 


C 


CDOI 


ACLOI = zero lift coefficient angle of attack 
ASIP = stall angle of attack 

CDBSI = drag coefficient at stall angle of attack smooth, 18% 
thickness to chord ratio and Reynolds number of 3 x 10^. 

CDOI = minimum drag coefficient 

RCL = reduction in lift coefficient at stall angle 

LEXP = power of the curvefit data after stall 

SLI = slope of lift coefficient 

ASCLOI = a constant, 90° 

RENS = Reynolds number 

FTL = slope correction factor for thickness to chord ratio and the 
Reynolds number effect. 

FTD = drag coefficient correction for thickness to chord ratio and 
Reynolds number. 

BCLF = zero incident angle of attack 

ASFP = stall angle of attack corrected for thickness to chord 
ratio and Reynolds number. 

ASFN = same as ASFP, but negative. 

CLSPF = maximum lift coefficient 
CLSNF = lift coefficient at ASFN 
CDMAX = maximum drag coefficient 

Secondary subroutines called: 

None 
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Notes : 


Airfoil characteristics are based on a two-dimensional, smooth, 18% 
thickness to chord ratio airfoil at a Reynolds number of 3 x 10^. This 
data IS corrected for sectional thickness to chord ratios and Reynolds 
numbers. 


CD44 

This subroutine is the same as CD230 except that it computes airfoil 
characteri sties of the NACA 44XX airfoil series. 

CTRWT 


Function: The purpose of this subroutine is computation of the negative 

torque of a counter-weight for a single bladed wind turbine. 

Arguments: All arguments are transmitted by COMMON except 

V = wi nd speed 
RHB = hub radius 

CTRQ = negative torque of counter-weight 


Secondary subroutines called: 
SPLINT 
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FACTORS 


Function: The purpose of this subroutine is the calculation of all 

induction factors. 

Arguments: All arguments are transmitted by COMMON except 

ZETA = Radial coordinates 
NP = number of points used in the spline fit 
J = index referring to the corresponding point on the blade 
XHB = dimensionless hub radius 
YI = induction factor 
YIBZ = induction factor at ZETP-DEL 
YIAZ = induction factor at ZETP+DEL 
SKPl) 

} = interval around the singularity 
SKP2) 

Secondary subroutines called: 

CALCI, SPLCOE 
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GAUSS 


Function: The purpose of this subroutine is the solution of a linear system 

of equations using the Gauss elimination technique. 


Arguments: All arguments are transmitted by COMMON except 

CA = coefficient matrix 
Ml = number of equations 
RQ = known vector 


Local variables: 

EPS = tolerance of the Gauss elimination technique 


Secondary Subroutines called: 
None 


Not es: 

This subroutine contains the following error message "error due to 
incorrect input of the number of equations". 
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GLQUD 


Function: The purpose of this subroutine is for setting up the Gauss- 

Laguerre integration. 

Arguments: All arguments are transmitted by COMMON except 

AUX = auxiliary subroutine containing the integrand 
ANS = computed integral 

Local Variables: 

X = Gauss-Laguerre roots 
W = weighting factor 

Secondary subroutines called: 

AUX 

SKIP 

Function: The purpose of this subroutine is the calculation of the interval 

around the singularity where the induction factors are 
1 nterpolated. 

Arguments: All arguments are transmitted by COMMON except 

ZETP = dimensionless radial position 
VEL = wind velocity 
OMEGA = rotational speed 
RB = blade radius 
SKPl ) 

} = interval around the singularity 
SKP2 ) 
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Secondary subroutines called: 

None 
Notes : 

The intervals around the singularities are a part of input data to the 
program. This subroutine is provided to facilitate the evaluation of 
proper intervals. The procedure, however, is more based on the 
numerical experimentation than theoretical analysis. Unrealistic 
results may occur if improper intervals are used, therefore, care must 
be taken in application of this subroutine for tip speed ratios outside 
the range of cases considered herein. For such cases, the proper 
intervals can be found from the fact that the induction factor, by 
definition, is a smooth and continuous function of radial position. 


SPLINT 

Function: The purpose of this subroutine is the computation of definite 

i ntegrals. 

Arguments: All arguments are transmitted by COMMON except 

NO = number of data points 
ZZ = spanwise location 
C = arrciy containing the integrand 
SUM = computed integral 

Secondary subroutines called: 

SPLCOE 
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SPLCOE 


Function: The purpose of this subroutine is the calculation of the 

coefficients for the spline fit. 

Arguments: All arguments are transmitted by COMMON except 

XP = vector containing independent variable 
YP = vector containing dependent variable 
M = number of data points 

Secondary subroutines called: 

None 


TINGRL 

Function: The purpose of this subroutine is the calculation of the 

integral part of the induced velocity. 


Arguments: All arguments are transmitted by COMMON except 

NS = index referring to the corresponding spanwise station 
NP = number of points on the blade for computing the integral 
ZETA = non-dimensional radial coordinate 
XHB = non-dimensional hub radius 
DELI = a small number, 10“^ 

YI = induction factor 

YIBZ = induction factor at ZETP-DEL 

YIAZ = induction factor at ZETP+DEL 

TINT = integral part of the induced velocity 

Secondary subroutines called: 


SPLINT 
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COMMON AND I/O PARAMETERS 


In this section, several important program details not covered in a 
description of all subprograms and a listing of program variables will be given. 
In particular, there will be a listing of terms that appear in COMMON followed 
by description of the input and output parameters of the program. 

COMMON 

Table 1 is a visual display of the subroutines in which COMMON appears. 


COMMON 

Subroutines where the common is used 

AAA 

AUX, CALI 

BBB 

MAIN, ASSGN, GLQUD 

CCC 

MAIN, ASSGN, CINTL, AUX, ALPHAE, TINGRL, 
CALPWR, CD230, CD44, FACTRS 

ODD 

MAIN, CD230, CD44, CTRW, CALPWR 

EEE 

MAIN, ASSGN, CINTL, ALPHAE, CD230, CD44 

FFF 

MAIN, ASSGN, CINTL, ALPHAE, CALPWR 

GGG 

MAIN, ASSGN, CALPWR 

HHH 

MAIN, ASSGN 

000 

ASSGN, CTRWT 


TABLE 1 
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COMMON AAA: 


contains the following parameters: 

KK = index of loop on the number of blades 

COMMON BBB: 

contains the following parameters: 

X = roots of Laguerre polynomial 
W = weight factors for the Gauss-Laguerre quadrature 

COMMON CCC: 

contains the following parameters: 

TT = radial location on the blade 
TIP = radial location on the blade 
PI = constant n 
EL = tip speed ratio 
NB = number of blades 

COMMON ODD: 

contains the following parameters: 

ROC = radian to degrees conversion 
WO = rotational speed 
RHO = density 

COMMON EEE: 

contains parameters that are described in the input data section which 

follows. They remain unchanged during the execution of the program. 
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COMMON FFF: 


Most parameters in this common are described in the following section on 
input parameters. The remaining parameters are: 

ZETP = dimensionless radial coordinate along the blade 
(equally spaced from hub to tip) 

BETA = twist angle distribution 

COMMON GGG: 

This common contains constants, calculated in the ASSGN subroutine, which 
are used throughout the main program and subroutines. These parameters are: 
CSSI = cosine of the coning angle 
COEF = a constant defined by 1/2 pRN 
COEFl = a constant defined by tt/ 2 pR2 


COMMON HHH: 

This common contains parameters that are described in the input data section. 


COMMON 000: 

This common contains information about the geometry of counter-weight of 
a single bladed rotor such as dimensions, coning angle, etc. All parameters 
are described in the section on input parameters. 
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DESCRIPTION OF THE INPUT PARAf€TERS 


lOP = output level control parameter 

0: output does not include the input data 

1: input data is first printed out then computed parameters are. 

ISKP = control parameter 

0: interpolation interval is input 

1: interpolation interval is calculated 

ITEM = control parameter 

1: wind speed variable 

2: rotational speed variable 

IREN = control parameter 

0: Reynolds' number effect not considered 

1: Reynolds' number effect considered 

MAXITR = maximum number of iterations 

NCASES = number of problem cases considered in the computer run 

N = number of stations along the blade. A maximum of 9 points is 

allowed (due to the machine storage limitation) 
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NPROF 


= parameter defining the type of airfoil being used. The present 
program is equipped with routines that can handle only two types 
of airfoils, NACA 230XX and NACA 44XX series. The following 
options are allowed: 

NPROF = 23000 (230XX series) 

NPROF = 4400 (44XX series) 

AO = slope of lift coefficient per degree. This parameter should be 
input as accurately as possible because of the dependence of the 
initial guess of the iteration loop upon it. 

BO = zero incident lift coefficient 

DELT = increment of the loop on NCASES 

NB = number of blades 

RB = radius of the blades 

OMEGA = rotational velocity, RPM. Initial rotational velocity if ITEM = 2. 
VEL = wind speed, mps. Initial wind speed if ITEM = 1. 

SI = coning angle, degree 

TCR75 = thickness to chord ratio at 3/4 of span. 
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SLTCR 

CI75 

XX 

CHRD 

BT 

X 

W 


slope of the thickness to chord ratio (TCR) distribution. This 
implies that approximately linear distribution of TCR can be 
handled. However, non-linear distributions can be easily handled by 
replacing the linear equation for TCR calculation with the 
non-1 1 near one. 

= chord at 3/4 of the span 

arrc^ containing the spanwise stations where information about 
chord and twist angle is provided. This input data must be provided 
in dimensionless radius. The first station must be located where 
blade begins i.e. , the dimensionless hub. 

array containing chord distribution. This data is to be provided 
at each XX location. 

twist angle distribution. This data is to be provided at each XX 
1 ocation. 

array containing roots of Laguerre polynomial. 

array containing weighting factor for Gauss-Laguerre integration. 
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Conditional Input Data 


The following parameters are input when ISKP = 0 

SKPl, SKP2 = interval around the singularity where the induction 
factors are interpolated. 

The following parameters are input when NB = 1 

(information about the size and location of different components of 
the counter-weight assuming an ellipsoidal counter-weight and 
cylindrical support span). 

RA = radius at which counterweight is located 

B1 = largest dimension of support spar 

B2 = smallest dimension of support spar 
B3 = dimension of counterweight (normal) 

B4 = dimension of counterweight (in radial direction) 

SIP = coning angle of spar 

N2 = number of points on support spar for integration 

N3 = number of points on counter-weight for integration 
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Description of the output. 


The output is composed of two parts: 

1) information pertaining to blade and operating condition. 
This includes 

- number of blades, blade radius, rpm 

- pitch and chord distribution 

- integration constants 

2) distribution of computed parameters along the blade such as 

- ci rculation 

- induced angle of attack 

- effective angle of attack 

- axial force 

- torque 

followed by the performance parameters 

- rotor power 

- alternator power 

- power coefficient 

- axial force coefficient 
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A SAMPLE PROBLEM 


The purpose of this section is to illustrate the use of VORTEX by way of an 
example problem. No attempt was made to select a problem that was trivial. 

Rather, the problem chosen was used because it is typical of most work done thus 
far. All input parameters needed to run VORTEX are displayed. The program output 
is presented as it comes from the machine. This output can serve as a means of 
deciding if the program is running correctly. 

INPUT 


The input data supplied the program is shown in Table 2. Explanation of this 
table is as follows: 


The first row of the table is the input data called from the first READ 
statement in the subroutine ASSGN. The parameters read are in order: 


lOP = 1 
ISKP = 1 

ITEM = 1 

IREN = 0 
MAXITR = 28 

NCASES = 1 


(the input data will be printed before the computed parameters 
(the interpolation interval will be determined within the 
subroutine SKIP) 

(the calculations will be performed for different wind speeds 
for a fixed rotational speed) 

(no Reynolds number correction of the airfoil data will be made 
(a maximum of 28 iterations are allowed to converge the 
circulation distribution) 

(only one problem is being run) 
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N =9 (there are 9 stations along the blade) 

NPROF = 23 (the blade is made of a NACA 230XX airfoil) 

DELT = 1.0 (if more than one problem were considered, the wind velocity 
[since ITEM = 1] would be increased by 1.0 mps) 

The second row of Table 2 is the input data called from the second READ 
statement in the subroutine ASSGN. The parameters read are in order: 


NB 

= 2 

(the machine considered has two blades) 


RB 

= 14.00 

(the 

blade radius in meters) 


OMEGA 

= 40.0 

(the 

rotor RPM) 


VEL 

= 12.0 

(the 

wind velocity in meters per second) 


SI 

= 3.00 

(the 

coning angle in degrees) 


AO 

= 0.0851 

(the 

slope of the lift coefficient curve 

per degree) 

BO 

= 0.1030 

(the lift coefficient at zero degrees of 

incidence) 


The third row of Table 2 is the input data called from the third READ 
statement in the subroutine ASSGN. The parameters read are in order: 


TCR75 = 0.240 

SLTCR = 0.0920 

CI75 = 1.520 


(the thickness to 
(the slope of the 
chord ratio) 

(the chord length 


chord ratio at 3/4 i 
distribution of the 

in meters at 3/4 of 


f the blade span 
thickness to 

the blade span) ' 
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1 1 

1 0 28 

2 14.00 

40.0 

.240 

.0920 1.520 

0.31648D0 

OO.OODO 

0.40008D0 

OO.OODO 

0.50000DO 

OO.OODO 

0.70000D0 

OO.OODO 

0.80000D0 

00.0000 

0.82912D0 

OO.OODO 

0.88608D0 

OO.OODO 

0.94304D0 

OO.OODO 

l.OOOOODO 

OO.OODO 


.81498279233948890D2 

.69962240035105030D2 

.61058531447218762D2 

.53608574544695070D2 

.47153106445156323D2 

.41451720484870767D2 

.36358405801651622D2 

.317760413523747-23D2 

.27633971743327170D2 

.23887329848169733D2 

.20491460082616425D2 

.17417992646508979D2 

.14642732289596674D2 

.12146102711729766D2 

.99120980150777060D1 

.79275392471721520D1 

.61815351187367654D1 

.46650837034671708D1 

.33707742642089977D1 

.22925620586321903D1 

.14255975908036131D1 

•76609690554593660D0 

.31123914619848373D0 

.59019852181507977D-1 


9 23 1.0 

3.00 0.0851 0.1030 

1.5200D0 
1.5200D0 
1.5200D0 
1.5200D0 
1.5200D0 
1.4230D0 
1.2700D0 
1.2200D0 
1.1700D0 
.55753457883283568D-34 
.40883015936806578D-29 
.24518188458784027D-25 
.36057658645529590D-22 
.20105174645555035D-19 
.53501888130100376D-17 
.78198003824594480D-15 
. 689418 10529580857D- 13 
.39177365150584514D-11 
. 15070082262925849D-09 
.40728589875499997D-08 
.79608129591336300D-07 
.11513158127372799D-05 
. 12544721977993333D-04 
. 10446121465927518D-03 
.67216256409354789D-03 
.33693490584783036D-02 
.13226019405120157D-01 
.40732478151408646D-01 
.98166272629918890D-01 
.18332268897777802D0 
.25880670727286980D0 
.25877410751742390D0 
.14281197333478185D0 


TABLE 2 

INPUT PARAMETERS 
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The fourth through twelfth rows in Table 2 are the input data called from 
the fourth READ statement in the subroutine ASSGN. The parameters read are in 
order: 

XX(J), J = (the first column; 9 spanwise locations along the 

blade: XX(9) = 1.000 is the blade tip) 

BT(J), J = 1,...,9 (the second column; 9 values of blade twist angle) 

CHRD(J), J = 1,...,9 (the third column; 9 values of the blade chord length 

in meters) 

The last 24 rows in Table 2 are the input data called from the fifth READ 
statement in the subroutine ASSGN. The parameters read are in order: 

X(L), L = 1,...,24 (the first 24 roots of the Laguerre polynomials) 

W(L), L = 1,...,24 (the first 24 weighting factors to be used in the 

Gauss-Laguerre integration) 
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Output 


Table 3 is a presentation of the output data of VORTEX. It can be seen that 
all input data has been written before the calculated parameters. After which the 
following arrc^ys are presented in order: 


nondimensional radial blade location 

= ZETP(IH) 

ci rculation 

= GAMM(IH) 

induced angle of attack 

= ALI(IH) 

effective angle of attack 

= ALE(IH) 

axial force on rotor 

= AFRCE(IH) 

rotor torque 

= TRQ(IH) 


The last row in Table 3 is a presentation of the computed output data. 


VEL = 

12.0 

(wind speed, meters/second) 

OMEGA = 

40.0 

(blade rotational speed, RPM) 

XRA = 

4.8802 

(tip speed ratio: fiR/V) 

SUMQ = 

42812.612 

(total torque on rotor. Newtons) 

RPWR = 

179.333 

(rotor power, kw) 

APWR = 

170.366 

(alternator power, kw) 

CP = 

0.276 

(power coefficient) 

SOME = 

26290.120 

(total axial force on rotor. Newtons) 

CF = 

0.485 

(axial force coefficient) 
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LOCATION 
0.31648E+00 
0.40192E+00 
0.48736E+00 
0.57280E+00 
0.65824E+00 
0.74368E+00 
0.82912E+00 
0.91456E+00 
0. 10000E+01 

MRS RPM 

12.0 40.0 


OPERATING CONDITIONS 


NUMBER OF BLADES 

2 

RADIUS OF BLADE, m 

14.00 

CONING ANGLE, degree 

3 0 

ROTATIONAL SPEED, rpm 

40 0 

WIND SPEED, m/s 

12.0 

TYPE OF AIRFOIL 

23 



BLADE DATA 


LIFT 

COEFF. SLOPE 

0.085 

ZERO 

INCIDENT LIFT 

0. 103 

NUMBER OF STATIONS 

9 

THICKNESS @ 3/4 SPAN 

0 24 

THICKNESS DIST. SLOPE 

0.09 

CHORD 

@ 3/4 SPAN 

1.52 

LOCATION 

TWI ST( degree) 

CHORD(m) 

0.316480 

0.000000 

1.520000 

0.400080 

0 000000 

1.520000 

0.500000 

0.000000 

1 520000 

0.700000 

0.000000 

1.520000 

0.800000 

0.000000 

1 520000 

0.829120 

0.000000 

1.423000 

0 886080 

0.000000 

1 .270000 

0.943040 

0.000000 

1.220000 

1.000000 

0.000000 

1.170000 


INTEGRATION CONSTANTS 


0.81498279233948890+02 
0.69962240035105030+02 
0.61058531447218760+02 
0.53608574544695070+02 
0.47153106445156320+02 
0.41451720484870770+02 
0.36358405801651620+02 
0.31776041352374720+02 
0.27635971743327170+02 
0.23887329848169730+02 
0 . 204914600826 1 6420+02 
0.17417992646508980+02 
0. 14642732289596670+02 
0. 12146102711729770+02 
0.9912098015077706D+01 
0.79275392471721520+01 
0.61815351187367650+01 
0.46650837034671710+01 
0 33707742642089980+01 
0.22925620586321900+01 
0. 14255975908036130+01 
0.76609690554593660+00 
0.31123914619848370+00 
0 . 5901985218 1 507980-0 1 


0.55753457883283570-34 
0 . 40883015936806580-29 
0.24518188458784030-25 
0.36057658645529590-22 
0.20105174645555030-19 
0.53501888130100380-17 
0.78198003824594480-15 
0.68941810529580860-13 
0.39177365150584510-11 
0. 15070082262925850-09 
0.40728589875500000-08 
0.79608129591336300-07 
0 11513158127372800-05 
0. 12544721977993330-04 
0 1044612146592752D-03 
0.67216256409354790-03 
0 33693490584783040-02 
0.13226019405120160-01 
0.40732478151408650-01 
0.98166272629918890-01 
0. 18332268897777800+00 
0.25880670727286980+00 
0.25877410751742390+00 
0.14281197333478180+00 


CIRCULATION 
0 OOOOOE+00 
0.21904E+02 
0 26060E+02 
0.33172E+02 
0 37414E+02 
0. 38242E+02 
0.37602E+02 
0.28642E+02 
-0. 13419E-14 

XRATIO 

4.8802 


ALPHA I 
-0.59498E+00 
-0. 13449E+00 
-0 18203E-01 
-0.93133E-01 
-0 52209E-01 
-0.82225E-01 
-0.67607E-01 
-0.83125E-01 
-0.22293E+00 

TORQUE 

42812.612 


ALPHA 0 
-0.21021E-01 
0 33643E+00 
0.37931E+00 
0.24998E+00 
0.24919E+00 
0. 18629E+00 
0. 17436E+00 
0. 13699E+00 
-0.21086E-01 


A. FORCE 
-0.67172E+01 
0. 18568E+05 
0.27341E+05 
0 36660E+05 
0 48136E+05 
0.53644E+05 
0.60018E+05 
0 53593E+05 
-0 67986E+02 


TORQUE 

-0.49850E+03 
0 14388E+05 
0.24108E+05 
0.63859E+05 
0.98668E+05 
0.95743E+05 
0. 1 1 143E+06 
0.83884E+05 
-0.79514E+04 


POWER ALT POWER CP A. FORCE CF 

179.333 170.366 0 276 26290.120 0 48 


TABLE 3 

OUTPUT PARAMETERS 
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CLOSURE 


In this report, a computer program entitled VORTEX for the aerodynamic 
performance prediction of horizontal axis wind turbines has been presented. 

The program is fairly general as it can handle wind turbines with any number of 
blades having a variety of blade geometry and airfoil section. It has been 
found that the program VORTEX is relatively efficient: most problems are 

solved in less than a minute of IBM 4341 CPU time. Predicted results have been 
found in good agreement with existing experimental data. 
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PROGRAM LISTING 


SUBROUTINE AAA 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

THIS PROGRAM CALCULATES THE THEORETICAL AERODYNAMIC PERFORMANCE 
PARAMETERS OF A HORIZONTAL AXIS WIND TURBINE USING A HELICAL 
VORTEX METHOD. THE CONTROL PARAMETERS ARE: 

IOP=l STANDARD OUTPUT 

IOP=0 STANDARD AS WELL AS INPUT DATA 

ITEM=1 WIND SPEED IS VARIABLE 

ITEM=2 ROTATIONAL SPEED IS VARIABLE 


THE IMPORTANT INPUT PARAMETERS ARE: 


AO 

BO 

DEL 

DELT 

EL 

EPS 

NCASES 

N 

NB 

NP 


OMEGA 

PR 

RB 

RA 

SI 

SIP 

SKP1( J) 
SKP2( J) 
VEL 
X,W 
ZETP 


LIFT CURVE SLOPE 

ZERO INCIDENT LIFT COEFFICIENT 

THE INTERVAL AROUND SINGULARITY THAT IS INTEGRATED 
SEPARATELY. 

THE INCREMENTAL CHANGE IN PARAMETRIC RUN 
THE TIP SPEED RATIO, V/WO*R 

A SMALL NUMBER, TOLERANCE OF GAUSS ELIMINATIO 
METHOD 

NUMBER OF CASES IN PARAMETRIC RUN 
NUMBER OF POINTS ON THE BLADE 
NUMBER OF BLADES 

NUMBER OF POINTS USED IN THE SPLINE INTERPOLATION 
AND INTEGRATION. IT MUST BE AT LEAST 4*N SO THAT 
WHEN INTEGRATING FROM X=0 TO ZETP-DEL OR FROM ZETP+DEL 
TO 1, THERE WILL ALWAYS BE AT LEAST 4 POINTS FOR 
THE SPLINE INTEGRATION. (SUBROUTINE SPLCOE WILL 
FAIL WITH LESS THAN FOUR POINTS) 

ROTATIONAL VELOCITY, RPM 
RATED POWER OF WIND TURBINE, KW 
ROTOR RADIUS, FT 

RADIUS AT WHICH COUNTER- WEIGHT IS LOCATED 
CONING ANGLE, DEGREES 

CONING ANGLE OF SPAR SUPPORT, DEGREES 
BLADEWISE LOCATIONS BETWEEN WHICH THE INDUCTION 
FACTOR IS INTERPOLATED 
WIND VELOCITY, MPH 

LAGUERRE-GAUSSIAN INTEGRATION CONSTANTS. 
NON-DIMENSIONAL DISTANCE ALONG THE BLADE WHERE 
CALCULATIONS HAVE BEEN PERFORMED FOR FOURIER 
SERIES COEFFICIENTS. 


MAIN PROGRAM 


IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION ALE(9), ALI(9), RAS(9), 

1 ALG(9 ),ZETA(81), 

2 CA(9 ,9), AS(9), TIGRL(9,9), 
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3 GAMM(9) ,AFRCE(9) ,TRQ(9) 

5 ,SKP1(9) ,SKP2(9) 

C0MM0N/BBB/X{24) , W(24) 

C0MM0N/CCC/P2 , TT , TTP , P I , EL , NB 

COMMON/DDD/ RDC,WO,RHO 

COMMON/EEE/ RB, TCR75 , SLTCR, CI75 

COMMON/FFF/ N , NPROF , AO , BO , 2ETP ( 9 ) , BETA ( 9 ) , C ( 9 ) 

COMMON/GGG/ CSS I , COEF , COEFl 
COMMON/HHH/ MAXITR,NCASES, ITEM, BELT 
DATA NCASES,DELT, ITEM/5, 2. DO, 1/ 

DATA IN, 10/3,8/ 

OPEN (NAME=* RTPM.DAT' ,TYPE='OLD’ , READONLY, UNI T= IN) 

RADIANS TO DEGREES CONVERSION 

READ IN AND PRINT OUT NECESSARY PARAMETERS 

CALL ASSGN ( OMEGA , VEL , P I , NB , R , SKP 1 , SKP2 , RDC ) 

NUMB=1 

NP=(N-1)*10+1 
DO 10 1=1, NP 

10 ZETA ( I ) =FLOAT ( I - 1 ) /FLOAT ( NP- 1 ) * ( 1 . - ZETP ( 1 ) ) + ZETP ( 1 ) 

20 V=VEL 

WO=OMEGA*P I /3 0 . DO 
EL=V/RB/WO 

CALCULATE THE INITIAL VALUE FOR A'S 
CALL Cl NTL ( SKP 1 , SKP2 , NP , ZETA , RDC , ALG , R , WO , T I GRL 
1,CA,RAS) 

THIS SECTION COMPUTES ALE VALUES USING AN ITERATION PROCESS 
CALL ALPHAE ( CA , RAS , T I GRL , R , WO , MAX I TR , ALG , AL I , ALE , GAMM ) 


CALL CALPWR ( ALE , R , RB , V , AFRCE , TRQ , SUMF , SUMQ , CF , CP , RPWR 
1,APWR,XRA) 


PRINT OUT STANDARD OUTPUT 

CALL OUTPUT ( N , ZETP , GAMM , ALE , AL I , AFRCE , TRQ , VE L , OMEGA , XRA , SUMQ 
1 , RPWR , APWR , CP , SUMF , CF ) 

NUMB=NUMB+1 

IF(NUMB.GT.NCASES) GO TO 30 
IF(ITEM.EQ.l) VEL=VEL+DELT 
I F ( I TEM . EQ . 2 ) OMEGA=OMEGA+DELT 
GO TO 20 
30 CONTINUE 

CLOSE ( UNIT=IN ) 

STOP 

END 
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SUBROUTINES OF THE PROGRAM 

SUBROUT I NE AS SGN ( OMEGA , VEL , P I , NB , R , SKP 1 , SKP2 , RDC ) 

THIS SUBROUTINE READS IN AND PRINTS OUT THE INPUT DATA 

THE OUTPUT PARAMETERS IN THE ARGUMENT ARE: 

OMEGA ROTATIONAL VELOCITY, RPM 

VEL WIND SPEED 

SI CONING ANGLE, DEGREES 

RDC RADIANS TO DEGREES CONVERSION 

MAXITR MAXIMUM NUMBER OF ITERATIONS 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION SKP1(9) ,SKP2(9) ,CHRD(9) ,BET(9) ,XX(9) 
COMMON/BBB/X(24) , W(24) 

COMMON/EEE/ RB, TCR75, SLTCR, CI75 
COMMON/FFF/ N , NPROF , AO , BO , ZETP ( 9 ) , BETA ( 9 ) , C ( 9 ) 

COMMON/GGG/ CSS I , COEF , COEFl 
COMMON/HHH/ MAXITR, NCASES, ITEM, DELT 
COMMON/000/ RA , B 1 , B2 , B3 , B4 , S I P , N2 , N3 
DATA IN, 10/3,8/ 

RDC=57 . 2957795130823D0 
RHO=1.2254710D0 
PI=3 . 14159265358979D0 

READ (IN, 100) IOP,ISKP,I TEM , I REN , MAX I TR , NCASES , N , NPROF , DE LT 
READ ( I N , 1 1 0 ) NB , RB , OMEGA , VE L , S I , AO , BO 
IF(IOP.EQ.l) WRITE (10,120) NB, RB, SI , OMEGA, VEL, NPROF 
READ (IN, 130) TCR75,SLTCR,CI75 

IF(IOP.EQ.l) WRITE(IO, 150) AO, BO, N, TCR75 , SLTCR, CI75 

READ IN BLADE ANGLE, CHORD, AND CORRESPONDING LOCATIONS 
ON THE BLADE 

DO 10 J=1,N 

READ(IN,160) XX(J), BT, CHRD(J) 

IF(IOP.EQ.l) WRITE(I0,170) XX( J) ,BT,CHRD( J) 

BET(J)=BT/RDC 
10 CONTINUE 

C READ IN DATA FOR A 24 POINT LANGUERRE- GAUSS INTEGRATION 

READ(IN,140) (X(L), W(L), L=l,24) 

IF(IOP.EQ.l) WRITE( 10,220) 

IF(IOP.EQ.l) WRITE(I0,230) (X( I ) , W( I ) , 1=1, 24) 

2ETP(1)=XX(1) 

C(1)=CHRD(1) 

BETA(1)=BET(1) 

J=2 

DO 20 1=2, N 

ZETP(I)=(1.D0-ZETP(1) )/FLOAT((N-l) )*(!-! )+ZETP(l) 

20 CONTINUE 
DO 40 1=2, N 

IF(ZETP(J) .LE.XX(I) .AND.ZETP(J) .GE.XX(I-l) ) GO TO 30 
GO TO 40 

30 DL=XX(I)-XX(I-1) 

DLP=ZETP( J)-XX(I ) 


51 



ooo ooo on 


C( J)=(CHRD( I )-CHRD( I-l ) )/DL*DLP+CHRD( I ) 

BETA ( J ) = ( BET ( I ) -BET ( I - 1 ) ) /DL*DLP+BET ( I ) 

J=J+1 

IF(ZETP(J) .LE.XX(I) .AND.ZETP(J) .GE.XX(I-l) ) GO TO 30 
40 CONTINUE 
62 DO 61 J=1,N 

C IF(IOP.EQ.l) WRITE(IO, 170) ZETP ( J) , BETA( J ) , C( J) 

61 CONTINUE 

IF(ISKP.EQ.l) GO TO 50 
READ(IN,12) (SKPl(J), SKP2 ( J ) , J=1 , N) 

GO TO 60 

50 CALL SKIP(ZETP, RB, OMEGA, VEL,N,SKP1,SKP2) 

60 CSSI=DCOS(SI/RDC) 

R=RB*CSSI 

COEF=0 . 5*NB*RB*RHO 
COEF1=0.5*RHO*PI*R**2 
AO=AO*RDC 

IF(NB.NE.l) GO TO 70 
READ IN COUNTER- WEIGHT DIMENSIONS 
SIP IS THE CONING ANGLE OF THE SPAR SUPPORT 
READ(IN,190) RA,B1,B2,B3,B4,SIP 

N2 A PARAMETER WHERE N2-1 IS THE NUMBER OF SUBDOMAINS IN SPAR 
SUPPORT FOR SPLINE INTEGRATION 

N3 SAME AS N2 EXCEPT THAT N3-1 IS FOR COUNTER- WEIGHT 
READ(IN,180) N2,N3 

IF(IOP.EQ.l) WRITE(IO,200) N2,RA, SIP,B1,B2 
IF(IOP.EQ.l) WRITE(IO,210) N3,B3,B4 
70 CONTINUE 

INPUT AND OUTPUT FORMATS 

100 FORMAT (8I5,F10.3) 

110 FORMAT(I5,6F10,5) 

120 F0RMAT(///41X, 'OPERATING CONDITIONS /41X, ' ' 

1,/35X, 'NUMBER OF BLADES' , lOX, 15, /35X, ' RADIUS OF BLADE, m' 
2,6X,F7.2,/35X, 'CONING ANGLE, degree ', 6X, F5 . 1 ,/35X, ' ROTATIONAL' 

3, ' SPEED, rpm' ,5X,F5.1,/35X, 'WIND SPEED, m/s' , 11X,F5. 1,/35X 

4, 'TYPE OF AIRFOIL' , IIX, 15///) 

130 FORMAT (3F10.4) 

140 FORMAT (2D30.17) 

150 FORMAT (45X, 'BLADE DATA',/45X, ' ',/35X, 'LIFT ' 

1, 'COEFF. SLOPE' ,9X,F5. 3, /35X, 'ZERO INCIDENT LIFT ' , 8X, F5 . 3 

2, /35X, 'NUMBER OF STATIONS ', 8X, 15 , /35X, ' THICKNESS ' 

3, ' @ 3/4 SPAN' ,6X,F5. 2, /35X, 'THICKNESS DIST. SLOPE' 
4,5X,F5,2,/35X, 'CHORD @ 3/4 SPAN' , lOX, F5 . 2 , //29X, ' LOCATION' 

5, 6X, ' TWI ST ( degree ) ' ,5X, 'CHORD(m) ' ) 

160 FORMAT{3F16.6) 

170 FORMAT(21X,3F16.6) 

12 FORMAT (2F10. 5) 

180 FORMAT (2 15) 

190 FORMAT(6F8.3) 

200 FORMAT ( /45X, ' SPAR DATA',/45X,' ' /35X, ' NUMBER OF STATIONS' 

1,7X, 15, /35X, 'RADIUS, m' , 25X, F5 . 2 , /35X, ' B1 DIMENSION, m',10X,F5.2, 
2/35X,'B2 DIMENSION, m',10X,F5.2) 
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210 FORMAT(/41X, 'COUNTER-WEIGHT DATA',/41X, ' ' 

1,/35X, 'NUMBER OF STATIONS' 7X, 15, /35X, ' B3 DIMENSION, m',10X,F5.2, 
2/3 5X, 'B4 DIMENSION, m',10X,F5.2) 

220 F0RMAT(/42X, ' INTEGRATION CONSTANTS' ,/42X, ' ' ) 

230 FORMAT(8X,2D35. 16) 

RETURN 

END 
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SUBROUT INK C I NTL ( SKP 1 , SKP2 , NP , ZETA , RDC , ALG , R , WO 
1,TIGRL,CA,RAS) 

THIS SUBROUTINE COMPUTE THE INITIAL VALUES OF (?A1 

THE INPUT PARAMETERS IN THE ARGUMENT ARE: 

NP NUMBER OF POINTS USED IN THE SPLINE INTERPOLATION 

ZETA RADIAL COORDINATES 

SI CONING ANGLE, DEGREES 

ALG GEOMETRIC ANGLE OF ATTACK 

R PARAMETER DEFINED BY RB*COS(SI) 

THE OUTPUT PARAMETERS IN THE ARGUMENT ARE: 

TIGRL MATRIX CONTAINING THE INTEGRAL FROM HUB TO TIP 
CA COEFFICIENT MATRIX <:A| 

RAS KNOWN VECTOR IN THE MATRIX EQUATION 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION IS(81) ,YI(81) ,TIGRL(9,9) ,CA(9,9) 

1,RAS(9) ,ALG(9) ,ZETA(81) 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 
COMMON/EEE/ RB, TCR75 , SLTCR, CI75 
COMMON/FFF/ N, NPROF , AO , BO , ZETP ( 9 ) , BETA ( 9 ) , C (9 ) 

DEL1=0.000001D0 
XHB=ZETP ( 1 ) 

DO 10 J=1,N 

ALG IS THE GEOMETRIC ANGLE OF ATTACK 
IF(ZETP( J) .LE.O.OOIDO) ALG( J)=-BETA( J) +PI/2 . ODO 
IF(ZETP( J) .GT.O.OOIDO) ALG( J)=-BETA( J) +DATAN(EL/ZETP( J) ) 

10 CONTINUE 
DO 40 J=1,N 
DO 20 1=1, NP 
IS(I)=0 

20 YI(I)=0.0D0 
TTP=ZETP( J) 

SH=EL* * 2 + ZETP ( J ) * * 2 

CALL FACTRS ( SKP 1 , SKP2 , I S , ZETA , NP , J , XHB , YI , YI BZ , YI AZ ) 

WRITE(7,2002) ( YI ( I ) , 1=1 ,NP ) 

2002 FORMAT (5F10. 4) 

DO 30 NS=1,N 

CALCULATE THE INTEGRAL FROM XHB TO 1 

CALL TINGRL ( NS , NP , ZETA , XHB , DELI , YI , YI BZ , YI AZ , T INT ) 

TIGRL( J,NS)=TINT 

ARG=P I * ( ZETP ( J ) -XHB ) / ( 1 . ODO-XHB ) 

T=- SLTCR* ( ZETP ( J ) -0 . 75 ) +TCR75 
SLC=AO+ ( 0 . 18-T ) *0 . 0050*RDC 
C SET UP THE MATRIX EQUATION 

CA( J,NS)=2.DO*R*DSIN(NS*ARG)/AO/C( J)-NS*C( J)*TINT/4.D0/(1.D0- 
1XHB) 

30 CONTINUE 
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T=-SLTCR* ( ZjETP ( J ) - 0 . 75 ) +TCR75 

RAS ( J ) =WO*R**2 *DSQRT ( SH ) * ( BO/AO+ALG( J ) + . 005* ( . 18-T ) /SLC ) *C ( J ) 
40 CONTINUE 
RETURN 
END 


55 



ooooonoooonon ooo 


SUBROUT I NE ALPHAE ( CA , RAS , T I GRL , R , WO , MAX I TR , ALG , AL I , ALE , GAMM ) 

THIS SUBROUTINE CALCULATES EFFECTIVE ANGLE OF ATTACK BY 

AN ITERATION SCHEME 

THE INPUT PARAMETERS IN THE ARGUMENT ARE: 

CA COEFFICIENT MATRIX (?A| 

RAS KNOWN VECTOR IN THE MATRIX EQUATION 

TIGRL MATRIX CONTAINING THE INTEGRAL FROM HUB TO TIP 
R PARAMETER DEFINED BY RB*COS(SI) 

MAXITR MAXIMUM NUMBER OF ITERATIONS 

ALG GEOMETRIC ANGLE OF ATTACK 

THE OUTPUT PARAMETERS IN THE ARGUMENT ARE: 

ALE EFFECTIVE ANGLE OF ATTACK 

ALI INDUCED ANGLE OF ATTACK 

GAMM CIRCULATION 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION AS(9) ,GAMM(9) , IS ( 81 ) , TIGRL( 9 , 9 ) , CA( 9 , 9 ) , RAS ( 9 ) 
1,ALG(9) ,ALE(9) ,ALI(9) 

COMMON/CCC/P2 , TT, TTP , PI , EL, NB 
COMMON/EEE/ RB, TCR75 , SLTCR, CI75 
COMMON/FFF/ N , NPROF , AO , BO , 2ETP ( 9 ) , BETA ( 9 ) , C ( 9 ) 

DATA 10/8/ 

RDC=57.2957795130823D0 
NITER=0 
XHB=ZETP(1) 

GO TO 40 
10 NITER=NITER+1 
DO 30 J=1,N 

VR=WO* R*DSQRT ( EL* *2 +ZETP ( J ) * *2 ) 

ARG=P I * ( ZETP ( J ) -XHB ) / ( 1 . ODO-XHB ) 

C REVISE MATRIX EQUATION FOR ITERATION PROCESS 

DO 20 NS=1,N 

IF(J.EQ.l .OR. J.EQ.N) CA( J, NS)=-TIGRL( J, NS ) *NS/4 . D0/( 1 . DO-XHB ) 
IF(J.EQ.l .OR. J.EQ.N) GO TO 20 
CA( J,NS)=DSIN(NS*ARG) 

20 CONTINUE 

IF(J.EQ.l .OR. J.EQ.N) GO TO 30 
RAS( J)=GAMM( J) 

30 CONTINUE 
40 DO 50 J=1,N 
50 AS(J)=RAS(J) 

C SOLVE MATRIX EQUATION USING GAUSS ELIMINATION TECHNIQUE 

CALL GAUSS (CA, N, AS) 

IF(NITER.NE.O) GO TO 80 
DO 70 NK=1,N 

VR=W0*R*DSQRT(EL**2+ZETP(NK) **2 ) 

ARG=P I * ( ZETP ( NK ) -XHB ) / ( 1 . ODO-XHB ) 

GAMM(NK)=O.ODO 
DO 60 K=1,N 

C CALCULATE THE CIRCULATION 

60 GAMM ( NK ) =GAMM ( NK ) + AS ( K ) *DS IN ( K* ARG ) 

T=-SLTCR* { ZETP ( NK ) -0 . 75 ) +TCR75 
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SLC=AO+ ( 0 . 18-T) *0 . 0050*RDC 
C CALCULATE THE EFFECTIVE ANGLE OF ATTACK 

ALE{NK)=2 . 0D0*GAMM(NK)/( AO*C(NK)*VR)-BO/SLC- . 005* ( . 18-T)/SLC 
LST=N+NK 

IF(ALE(NK) .LE.0.04100D0) IS(LST)=2 
IF(NK.EQ.l.OR.NK.EQ.N) GO TO 70 
IF(ALE(NK) .GE.0.25) ALE(NK)=ALG(NK)- . 01 
70 CONTINUE 
80 DO 140 NK=1,N 
LST=N+NK 
IS(NK)=0 
STIGR=0.0D0 

VR=WO*R*DSQRT(EL**2+ZETP(NK)**2) 

DO 90 K=1,N 

90 STIGR=STIGR+TIGRL(NK,K)*AS(K)*K 

IF(NITER.GE. 1 . AND . (NK. EQ . 1 . OR.NK. EQ. N) ) GO TO 100 
ALI(NK)=STIGR/(4.D0*R*VR*(1.D0-XHB) ) 

ALE ( NK ) =ALG ( NK ) + AL I ( NK ) 

100 AL=ALE(NK) 

XB=ZETP(NK) 

IF(NPROF .EQ. 44) GO TO 110 
CALL CD230( IREN,AL,XB,VR,CL,CD) 

GO TO 120 

110 CALL CD44(AL,CD,CL,XB,VR) 

120 CLO=2.DO*GAMM(NK)/VR/C(NK) 

DELCL= ( CL-CLO ) / ( CL+ . 00001 ) 

IF(DABS(DELCL) .GT. .02)GO TO 130 
GO TO 140 

130 IF(IS(LST) .EQ.2) GO TO 140 

IF(NK.EQ. 1. OR.NK. EQ.N) GO TO 140 
GAMM(NK)=.25DO*C(NK)*(CL+CLO)*VR 
IS(NK)=1 
140 CONTINUE 
SIA=0.D0 
DO 150 1=1, N 
150 SIA=SIA+IS(I) 

IF(SIA.LT.l) GO TO 170 
IF(NITER.GT.MAXITR) GO TO 160 
GO TO 10 

C ITERATIOIN PROCESS IS NOW COMPLETED 

160 WRITE( 10,200) 

STOP 

170 RETURN 

200 F0RMAT('=== NO. OF ITERATION EXEEDED THE LIMIT ===' ) 

END 
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c 

SUBROUTINE CALPWR ( ALE , R , RB , V , AFRCE , TRQ , SUMF , SUMQ , CF , CP 
1,RPWR,APWR,XRA) 

C THIS SUBROUTINE COMPUTES FORCE, TORQUE AND POWER 

C 

C THE INPUT PARAMETERS IN THE ARGUMENT ARE: 


c 

ALE 

EFFECTIVE ANGLE OF ATTACK 

c 

R 

PARAMETER DEFINED BY RB*COS(SI) 

c 

SI 

CONING ANGLE, DEGREES 

c 

RB 

RADIUS OF THE BLADE 

c 

V 

WIND SPEED 

c 

THE OUTPUT PARAMETERS IN THE ARGUMENT ARE 

c 

AFRCE 

LOCAL AXIAL FORCE ON THE BLADE 

c 

TRQ 

LOCAL TORQUE ON THE BLADE 

c 

SUMF 

TOTAL AXIAL FORCE 

c 

SUMQ 

TOTAL TORQUE 

c 

CF 

AXIAL FORCE COEFFICIENT 

c 

CP 

COEFFICIENT OF PERFORMANCE 

c 

RPWR 

TOTAL POWER 

c 

APWR 

ALTERNATOR POWER 

c 

XRA 

INVERSE OF TIP SPEED RATIO 


IMPLICIT 

REAL*8 (A-H,0-Z) 


DIMENSION ALE(9) ,AFRCE(9) ,TRQ(9) 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 
COMMON/DDD/ RDC,WO,RHO 

COMMON/FFF/ N, NPROF, AO, BO, ZETP( 9 ) , BETA( 9 ) , C( 9 ) 

COMMON/GGG/ CSSI , COEF, COEFl 
1600 DO 1750 NK=1,N 

PHI =BETA ( NK ) +ALE ( NK ) 

XB=ZETP(NK) 

VR=W0*R*DSQRT(EL**2+ZETP(NK)**2) 

ALOO=ALE(NK) 

IF(NPROF .EQ. 44) GO TO 1650 
CALL CD2 30(1 REN , ALOO , XB , VR , CL , CD ) 

GO TO 1700 

1650 CALL CD44(AL00,CD,CL,XB,VR) 

1700 CY=CL*DCOS(PHI ) +CD*DSIN(PHI ) 

CX=CL*DSIN(PHI ) -CD*DCOS(PHI ) 

COQT=COEF*C ( NK ) * VR* VR 

C THIS SECTION CALCULATES AXIAL FORCE AND TORQUE AT N LOCATIONS 

C ALONG THE BLADE 

TRQ ( NK ) =COQT*R* ZETP ( NK ) *CX 
1750 AFRCE(NK)=COQT*CY*CSSI 

C INTEGRATE THE AXIAL FORCE USING SPLINE INTERPOLATION FORMULA 

SUMF=0 . ODO 

CALL SPLINT(N, ZETP, AFRCE, SUMF) 

C INTEGRATE THE TORQUE OVER THE BLADE USING SPLINE INTERPOLATION 

C INTERPOLATION FORMULA 

SUMQ=0 . ODO 

CALL SPL I NT (N, ZETP, TRQ, SUMQ) 

C 

C CALCULATE THE EFFECT OF COUNTER- WEIGHT 
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c 

IF(NB.NE.l) GO TO 1850 
RHB=RB*ZETP(1) 

CALL CTRWT(V,RHB,CTRQ) 
1850 CF=SUMF/(C0EF1*V**2) 
SUMQ=SUMQ+CTRQ 
RPWR=SUMQ*WO/1000 . DO 
C COMPUTE ALTERNATOR POWER 

APWR= . 95D0* ( RPWR- . 07 5* PR) 
CP=SUMQ*W0/(C0EF1*V**3 ) 
XRA=WO*R/V 
RETURN 
END 
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SUBROUTINE OUTPUT ( N , ZETP , GAMM , ALE , AL I , AFRCE , TRQ , VEL , OMEGA , XRA , SUMQ 
1 , RPWR , APWR , CP , SUMF , CF ) 

THIS SUBROUTINE PRINTS OUT THE STANDARD OUTPUT 

THE INPUT PARAMETERS IN THE ARGUMENT ARE: 

N NUMBER OF TERMS IN THE SINE SERIES 

ZETP DIMENSIONLESS DISTANCE ALONG THE BLADE WHERE CALCULATION 

HAVE BEEN PERFORMED FOR FOURIER SERIES COEFFICIENTS 
GAMM CIRCULATION 

ALE EFFECTIVE ANGLE OF ATTACK 

ALI INDUCED ANGLE OF ATTACK 

AFRCE LOCAL AXIAL FORCE ON THE BLADE 
TRQ LOCAL TORQUE ON THE BLADE 

VEL WIND SPEED 

OMEGA ROTATIONAL VELOCITY, RPM 

XRA INVERSE OF TIP SPEED RATIO 

SUMQ TOTAL TORQUE 

RPWR TOTAL POWER 

APWR ALTERNATOR POWER 

CP COEFFICIENT OF PERFORMANCE 

SUMF TOTAL AXIAL FORCE 

CF AXIAL FORCE COEFFICIENT 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION ZETP(9) ,GAMM(9) ,ALE(9) ,ALI (9) ,AFRCE(9) ,TRQ(9) 

DATA 10/8/ 

WRITE(IO,2120) 

DO 1800 IH=1,N 

WRITE( 10,2130) ZETP(IH), GAMM(IH), ALI(IH), ALE(IH), AFRCE(IH), 

1 TRQ(IH) 

WRITE(IO,2140) 

WR I TE ( 1 0 , 2 1 5 0 ) VEL, OMEGA , XRA , SUMQ , RPWR , APWR , CP , SUMF , CF 

STANDARD OUTPUT FORMATS 

2120 FORMAT (//’ ',7X, ' LOCATION' , 8X, ' CIRCULATION' , 8X, ' ALPHA l' 

1 , lOX, 'ALPHA O', IIX, 'A. FORCE', IIX, 'TORQUE') 

2130 FORMATC ',6E17.5) 

2140 FORMAT (/7X, ' MPS ', 6X, ' RPM' , 9X, ' XRATIO ', 9X, ' TORQUE ', 12X, 

I'POWER* ,9X, 'ALT. POWER' ,5X, ' CP ' ,4X, 'A. FORCE' , 7X, 'CF* ) 

2150 FORMAT(2F10.1,F15.4,3F16.3,3F10.3,//) 

RETURN 

END 
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SUBROUT I NE CALC I ( X I , NOPT ) 

THIS SUBROUTINE CALCULATES INDUCTION FACTORS AS FOLLOWS: 
NOPT = 1 CALCULATES 11+12=1 

NOPT = 2 CALCULATES THE INTEGRAL OF M FROM 0 TO INFINITY 
WHERE I2=(ZETA-ZETP)*P2 
IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/AAA/ KK 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 
DIMENSION AN(2) 

EXTERNAL AUX 
IF(NB.EQ.l) GO TO 50 
DO 10 KK=1,2 

IF(NOPT.EQ.2 .AND. KK.EQ.2) GO TO 20 
CALL GLQUD(AUX,ANS) 

10 AN(KK)=ANS 

IF(NOPT.EQ.l .AND. DABS ( TT-TTP ) . LE . 0 . OOIDO ) XI=1.0D0 
IF(NOPT.EQ.l .AND. DABS ( TT-TTP ). GT. 0 . OOIDO ) 

1 XI=(AN(1)+AN(2) )*(TT-TTP) 

20 IF(NOPT.EQ.2) XI=AN(1) 

GO TO 100 

50 CALL GLQUD(AUX,ANS) 

ANK=ANS 

IF(NOPT.EQ.l .AND. DABS (TT-TTP ). LE . 0 . OOIDO ) XI=1.0D0 
IF(NOPT.EQ.l .AND. DABS (TT-TTP ). GT. 0 . OOIDO ) 

1 XI=ANK* (TT-TTP) 

100 RETURN 
END 
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SUBROUTINE GLQUD ( AUX , ANS ) 

C THIS SUBROUTINE SETS UP THE LAGUERRE- GAUSS INTEGRATION 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/BBB/X ( 2 4 ) , W ( 2 4 ) 

ANS=0 . ODO 
S=O.ODO 
DO 10 1=1,24 
Y=X(I) 

CALL AUX(Y,Z) 

10 S=S+Z*W(I) 

ANS=S 

RETURN 

END 
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SUBROUT I NE AUX ( THET , FX ) 

THIS SUBROUTINE CALCULATES THE INTEGRAND OF P1,P2&P3 
THE INTEGRAND OF P IS GIVEN BY A/B BELOW 
IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 
COMMON/AAA/ KK 
IF(NB.EQ.l) GO TO 10 
THETK=THET+2 . 0D0*PI* (NB-KK)/NB 

A=TT*TTP* ( TT-TTP*DCOS ( THETK ) ) + ( TT* ( THET*DS IN ( THETK ) +DCOS ( THETK ) ) 

1 -TTP)*EL**2 

B=(DSQRT(TT**2+TTP**2-2.0D0*TT*TTP*DCOS(THETK)+THET**2*EL**2) )**3 
1 *DSQRT(EL**2+TTP**2) 

GO TO 20 

10 A=TT*TTP* (TT-TTP*DCOS ( THET ) ) + ( TT* ( THET*DS IN( THET ) +DCOS ( THET ) ) -TTP ) 
1 *EL**2 

B=(DSQRT(TT**2+TTP**2-2 .ODO*TT*TTP*DCOS(THET)+THET**2*EL**2) )**3 
1 *DSQRT(EL**2+TTP**2) 

20 C=A/B 

FX=DEXP(THET)*C 

RETURN 

END 
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SUBROUTINE SPLCOE(XP, YP,M, CC) 

THIS SUBROUTINE CALCULATES THE COEFFICIENTS FOR SPLINE 

INTERPOLATION AS WELL AS INTEGRATION 

THE INPUT PARAMETERS IN THE ARGUMENT ARE: 

XP VECTOR CONTAINING INDEPENDENT VARIABLE 

YP VECTOR CONTAINING DEPENDENT VARIABLE 

M NUMBER OF DATA POINTS 

THE OUTPUT PARAMETER IN THE ARGUMENT IS: 

CC COEFFICIENT OF SPLINES 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION XP(81) ,YP(81) ,D(81),P(81),E(81),CC(4,81) , 
1AE(81,3) ,B(81) ,Z(81) 

MM=M-1 

DO 10 1=1, MM 
D(I)=XP(I+1)-XP(I) 

P(I)=D(I)/6. 

10 E(I)=(YP(I+1)-YP(I) )/D(I) 

DO 20 1=2, MM 
20 B(I)=E(I)-E(I-1) 

AE(1,1)=1.0D0 

AE(1,2)=-1.0D0-D(1)/D(2) 

AE(1,3)=D(1)/D(2) 

AE(2,3)=P(2)-P(1)*AE(1,3) 

AE(2,2)=2.0D0*(P(1)+P(2) )-P(l)*AE(l,2) 
AE(2,3)=AE(2,3)/AE(2,2) 

B(2)=B(2)/AE(2,2) 

DO 30 1=3, MM 

AE(I,2)=2.*(P(I-1)+P(I) )-P(I-l)*AE(I-l,3) 
B(I)=B(I)-P(I-1)*B(I-1) 

AE(I,3)=P(I)/AE(I,2) 

30 B(I)=B(I)/AE(I,2) 

QR=D(M-2)/D(M-l) 

AE(M, l)=1.0D0+QR+AE(M-2,3) 

AE(M,2)=-QR-AE(M, 1)*AE(M-1,3) 

B(M)=B(M-2)-AE(M, 1)*B(M-1) 

Z(M)=B(M)/AE(M,2) 

MN=M-2 

DO 40 1=1, MN 
K=M-I 

40 Z(K)=B(K)-AE(K,3)*Z(K+1) 

Z(1)=-AE(1,2)*Z(2)-AE(1,3)*Z(3) 

DO 50 K=1,MM 
QR=1.0D0/(6,0D0*D(K) ) 

CC(1,K)=Z(K)*QR 

CC(2,K)=Z(K+1)*QR 

CC(3,K)=YP(K)/D(K)-Z(K)*P(K) 

50 CC(4,K)=YP(K+1)/D(K)-Z(K+1)*P(K) 

RETURN 

END 
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SUBROUTINE GAUSS (CA, Ml , RQ) 

THIS SUBROUTINE SOLVES A SYSTEM OF EQUATIONS 

BY MEANS OF GAUSS ELIMINATION METHOD 

THE INPUT PARAMETERS IN THE ARGUMENT ARE: 

CA COEFFICIENT MATRIX <?Al 

Ml NUMBER OF EQUATIONS 

RQ KNOWN VECTOR IN THE MATRIX EQUATION 

THE OUTPUT PARAMETERS IN THE ARGUMENT ARE: 

RQ SOLUTION VECTOR 

lER ERROR MESSAGE 

EPS TOLERANCE OF THE GAUSS ELIMINATION TECHNIQUE 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION AO(81),RQ(9 ),CA(9,9) 

DATA EPS/0.5D0/ 

IEP=1 

DO 10 IP=1,M1 
DO 10 JP=1,M1 
AQ(IEP)=CA( JP, IP) 

IEP=IEP+1 
10 CONTINUE 

IF(M1)240,240,20 
20 IER=0 

PIV=0.0D0 

M3=M1*M1 

M4=M1 

DO 40 L=1,M3 
TB=DABS(AQ(L) ) 

IF(TB-PIV)40,40,30 
30 PIV=TB 
I=L 

40 CONTINUE 
TOL=EPS*PIV 

LST=1 

DO 180 K=1,M1 
IF(PIV)240,240,50 
50 IF(IER)80,60,80 
60 IF(PIV-TOL)70,70,80 
70 IER=K-1 
80 PIVI=1.0D0/AQ(I) 

J=(I-1)/M1 

I=I-J*M1-K 

J=J+1-K 

DO 90 L=K,M4,M1 
LL=L+I 

TB=PIVI*RQ(LL) 

RQ(LL)=RQ(L) 

90 RQ(L)=TB 

IF(K-M1)100,190,190 
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100 LEND=LST+M1-K 

IF(J)130,130,110 
110 II=J*M1 

DO 120 L=LST,LEND 
TB=AQ(L) 

LL=L+II 

AQ(L)=AQ(LL) 

120 AQ(LL)=TB 

130 DO 140 L=LST,M3,M1 
LL=L+I 

TB=PIVI*AQ(LL) 

AQ(LL)=AQ(L) 

140 AQ(L)=TB 

AQ(LST)=J 

PIV=0.0D0 

LST=LST+1 

J=0 

DO 170 II=LST,LEND 
PIVI=-AQ(II) 

IST=II+M1 

J=J+1 

DO 160 L=IST,M3,M1 
LL=L-J 

AQ(L)=AQ(L)+PIVI*AQ(LL) 
TB=DABS(AQ(L) ) 
IF(TB-PIV)160,160,150 
150 PIV=TB 
I=L 

160 CONTINUE 

DO 170 L=K,M4,M1 
LL=L+J 

170 RQ(LL)=RQ(LL)+PIVI*RQ(L) 
180 LST=LST+M1 


190 IF(M1-1)240,230,200 
200 IST=M3+M1 
LST=M1+1 
DO 220 1=2, Ml 
II=LST-I 
IST=IST-LST 
L=IST-M1 
L=AQ(L)+0.5D0 
DO 220 J=II,M4,M1 
TB=RQ(J) 

LL=J 

DO 210 K=IST,M3,M1 
LL=LL+1 

210 TB=TB-AQ(K)*RQ(LL) 
K=J+L 



RQ(J)=RQ(K) 

220 RQ(K)=TB 
230 GO TO 250 

240 WRITE(IO,300) 

250 RETURN 

300 F0RMAT(1 OX, 'ERROR DUE TO INCORRECT NUMBER OF EQUATIONS ( 
1 ZERO OR NEGATIVE'/) 

END 
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SUBROUT I NE CD2 3 0 ( I REN , ALPHA , X , W , CL , CD ) 

THIS SUBROUTINE CONTAIN THE LIFT AND DRAG AIRFOIL DATA 
CURVFFITS FOR NACA 230XXAIRFOIL( SMOOTH) WITH CORRECTIONS 
FOR THE EFFECT OF THICKNESS TO CHORD RATIO 

INPUT PARAMETERS ARE: 

IREN CONTROL PARAMETER 

ALPHA ANGLE OF ATTACK 

X SPANWISE LOCATION 

W RELATIVE VELOCITY 

OUTPUT PARAMETERS ARE: 

CL LIFT COEFFICIENT 

CD DRAG COEFFICIENT 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 
COMMON/EEE/RB, TCR75 , SLTCR, CI75 
COMMON/DDD/ RDC,WO,RHO 

DATA ACLOI , AS IP , CDBS I , CDOI/- 1 . 2D0 , 15 . DO , 0 . 02 15D0 , 0 . 0074D0/ 
DATA RCL,SLI,LEXP,ASCLO/.25, .103,5,90.0/ 

RENS=3000000.0 
A=ALPHA*RDC 
AALPHA=DABS ( ALPHA ) 

CALCULATE THE EFFECT OF THICKNESS TO CHORD RATO 


FR=X 

T=- SLTCR* ( FR- . 75 ) +TCR75 
IF(IREN.EQ.l) RENS=1000*W*CI75 
FTL= . 9534+2 . 03 ll*T-9 . 843700*T**2 

FTD=(1.0+2.0*(T-0. 18) )*(l-0.009*(0.000001*RENS-3.0) ) 
AR=RB/CI75 

CALCULATE LIFT CURVEFIT CONSTANTS 

IF(T.EQ.PRET) GO TO 20 

SLF=SLI*FTL 

BCLF=0 . -SLF*ACLOI 

IF (RENS.GT. 6000000. ) RENS=6000000 . 

ASFP=ASIP*(1+1. 1111*(0. 18-T) )*(l-0.02083*(3-.000001*RENS) ) 
ASFN=ACLOI-ASFP 

CLSPF=BCLF+SLF*ASFP-RCL* ( -55 . 3332*T**2+16 . 5332*T-0 . 1088 ) 
l+0.0275*(0.000001*RENS-3. ) 

CLSNF=BCLF+SLF*ASFN+RCL* ( -55 . 3332*T**2+16 . 5332*T-0 . 1088) 
l+0.0275*(0.0000001*RENS-3. ) 

SLSP=CLSPF/ ( ASFP- ASCLO ) 

SLSN=CLSNF/ ( ASFN+ASCLO ) 

BCLPS=CLSPF-SLSP*ASFP 

BCLNS=CLSNF-SLSN*ASFN 
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CALCULATE LIFT COEFFICIENT 

20 CONTINUE 

IF(A.LT.ASFN) GO TO 40 
IF(A.LT.O. ) GO TO 45 

IF(T.LT.0.18.AND.A.LT.ASFP*(l+0.02*(0.000001*RENS-3))) GO TO 501 

IF(T.LT.0.24.AND.A.LT.ASFP) GO TO 601 

IF(T.GE.0.24.AND.A.LT,ASFP) GO TO 701 

IF(T.LT.0.18.AND.A.LT.20. ) GO TO 502 

IF(T.LT.0.24.AND.A.LT.20. ) GO TO 602 

IF(T.GE.0,24.AND.A.LT.20. ) GO TO 702 

IF(A.LT.20. ) GO TO 55 

IF(A.LT.24. ) GO TO 57 

IF(A.LT.46) GO TO 59 

GO TO 60 

:0 CL=SLSN*A+BCLNS 
GO TO 80 

45 CL=SLF*A+BCLF+RCL* ( -55 . 3*T**2+16 . 5*T-0 . 109 ) * ( ABS ( A/ASFN) ) **LEXP 
l+(0.0275*(0.000001*RENS-3. ))* (ABS (A/ASFN) ) **LEXP 
GO TO 80 

501 CL=SLF*A+BCLF-RCL* ( -55 . 3*T**2+16 . 5*T-0 . 109+0 . 018* ( 0 . 000001* 
lRENS-3) )*(A/ASFP/(l+0.02*(0.000001*RENS-3) ) )**LEXP 
GO TO 80 

601 CL=SLF*A+BCLF-RCL* ( -55 . 3*T**2+16 . 5*T-0 . 109 ) * ( A/ASFP)**LEXP 
l+(0.0425*(0.000001*RENS-3. ) )* (A/ASFP)**LEXP 
GO TO 80 

701 CL=SLF*A+BCLF-RCL* ( -55 . 3*T**2+16 . 5*T-0 . 109 ) * ( A/ASFP ) **LEXP 
l+(0.01*(0.000001*RENS-3. ) )*(A/ASFP)**LEXP 
GO TO 80 

502 CL=(CLSPF-1 . 0 ) * ( A-20 )/( ASFP* ( 1+0 . 02* (0 . 000001*RENS-3 ) ) -20 ) +1 
GO TO 80 

602 CL=(CLSPF-1 . 0 ) * (A-20)/(ASFP-20) +1 . 0 
GO TO 80 

702 CL=(CLSPF-1 . 0-0 . 016667* ( 0 . 000001*RENS-3 ) ) * ( A-20 )/( ASFP-20 ) +1 . 0 
GO TO 80 

60 CL=1.5557*SIN(2*ALPHA) 

GO TO 80 
80 CONTINUE 

CALCULATE DRAG CURVEFIT CONSTANTS 

IF(T.EQ.PRET) GO TO 90 
CDBSF=CDBSI 

DF=(CDBSF-CDOI )/( ASFP-ACLOI ) **2 
CDMAX=1 . 5 

CDMAX=1 . 11+AR* . 0178 

DS=(CDMAX-CDBSF)/( 1 . 57-ASFP/RDC) **2 

CALCULATE DRAG COEFFICIENT 

90 IF(A.LT. (O.-ASFP) ) GO TO 100 

IF(T.GE.0.18.AND.A.LT.ASFP) GO TO 110 

IF(T.LT. .18.AND.A.LE.ASFP*(l+0.02*(0.000001*RENS-3) ) ) GO TO 111 
100 CD=CDMAX-DS*(1.57-AALPHA)**2 
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GO TO 130 

110 CD=CD0I*FTD+DF*(A-ACL0I)**2 
GO TO 130 

111 CD=CDOI*FTD+(CDBSF-CDOI ) * ( A-ACLOI ) **2/( ASFP* ( 1+0 . 02* ( 0 . 000001* 
lRENS-3) )-ACLOI)**2 

GO TO 130 
130 CONTINUE 
RETURN 
END 
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SUBROUTINE CD44 (ALPHA, CD, CL, X, W) 

THIS SUBROUTINE CONTAINS THE LIFT AND DRAG AIRFOIL DATA 
CURVFITS FOR NACA 44XX AIRFOIL( SMOOTH) WITH CORRECTIONS 
FOR THE EFFECT OF THICKNESS TO CHORD RATIO 

VARIABLE DEFINITIONS 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 
COMMON/DDD/ RDC,WO,RHO 
COMMON/EEE/ RB , TCR7 5 , S LTCR , C 1 7 5 

DATA ACLOI,ASIP,CDBSI,CDOI,PRET/-4.DO, 14. DO, .021D0, .0077D0, O.DO/ 
DATA RCL,SLI,ASCLO,LEXP/.33DO, . 0971D0 , 90 . DO , 5/ 

A=ALPHA*RDC 
AALPHA=DABS ( ALPHA ) 

CALCULATE THE EFFECT OF THICKNESS TO CHORD RATIO 


FR=X 

T=- SLTCR* ( FR- . 75 ) +TCR75 
AR=RB/CI75 

. . . CALCULATE THE EFFECT OF THICKNESS TO CHORD RATIO . . . 

FTL=-1 . 1329*T+1 . 2039 
FTD=1. 

. . . CALCULATE LIFT CURVFIT CONSTANTS . . . 

IF(T.EQ.PRET) GO TO 20 
SLF=1./(1./(SLI*FTL)+18.24/AR) 

BCLI=0 . - ( SLI *FTL ) *ACLOI 
BCLF=0 . -SLF*ACLOI 
CLSPI=(SLI*FTL)*ASIP+BCLI-SCL 
ASFP=ASIP*(1.+1.1905*( ,18-T)) 

ASFN=ACLOI-ASFP 

CLSPF=BCLF+SLF*ASFP-RCL* ( 5 . 5555*T**2-6 . 8433*T+1 . 8027 ) 
CLSNF=BCLF+SLF*ASFN+RCL* ( 5 . 5555*T**2-6 . 8433*T+1 . 8027 ) 
SLSP=CLSPF/( ASFP-ASCLO ) 

SLSN=CLSNF/ ( ASFN+ ASCLO ) 

BCLPS=CLSPF-SLSP*ASFP 
BCLNS=CLSNF- SLSN* ASFN 

. . . CALCULATE LIFT COEFFICIENT , . . 

20 CONTINUE 

IF(A.LT.ASFN) GO TO 40 
IF(A.LT.O. ) GO TO 45 
IF(A.LT.ASFP) GO TO 50 
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GO TO 60 

40 CL=SLSN*A+BCLNS 
GO TO 80 

45 CL=SLF*A+BCLF+RCL* ( 5 . 556*T*T-6 . 8433*T+1 . 8027 ) * (DABS( A/ASFN) ) **LEXP 
GO TO 80 

50 CL=SLF*A+BCLF-RCL* ( 5 . 556*T*T-6 . 8433*T+1 . 8027 ) * (DABS( A/ASFP ) ) **LEXP 
C 50 CL=SLF*A+BCLF-RCL*( A/ASFP )**LEXP 
GO TO 80 

60 CL=SLSP*A+BCLPS 
GO TO 80 
80 CONTINUE 

. . . CALCULATE DRAG CURVFIT CONSTANTS . . . 

IF(T.EQ.PRET) GO TO 90 
CDBSF=CDBSI 

DF=(CDBSF-CDOI )/ASFP**2 
CDMAX=1.11+.018*AR 
IF(AR.GT.50) CDMAX=2. 

DS=(CDMAX-CDBSF)/( 1 . 57-ASFP/RDC ) **2 

. . . CALCULATE DRAG COEFFICIENT . . . 

90 IF(A.LT. (O.-ASFP) ) GO TO 100 
IF(A.LT.ASFP) GO TO 110 
100 CD=CDMAX-DS*(1.57-AALPHA)**2 
GO TO 130 

110 CD=CD0I*FTD+DF*A**2 
GO TO 130 
130 CONTINUE 
PERT=T 

RETURN 
END 
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SUBROUT I NE CTRWT ( V , RHB , CTRQ ) 

C THIS SUBROUTINE CALCULATES THE NEGATIVE TORQUE OF COUNTER- WEIGHT 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/DDD/ RDC,WO,RHO 
COMMON/000/ RA , B1 , B2 , B3 , B4 , S I P , N2 , N3 
DIMENSION X2(81), Q(81), YI(81) 

SIPP=SIP/RDC 

DIP=DCOS(SIPP) 

Rl=RA-B4/2 . 0 

NP2=(N2-1)*3+1 

NP3=(N3-1)*5+1 

VEF=VEF*DIP 

ND=0 

DO 800 1=1, NP2 
ND=ND+1 

XZ(I)=(I-1)*(R1-RHB)/FL0AT(NP2-1) 

YI ( I )=B1- (B1-B2 ) *XZ( I )/(Rl-RHB) 

REFF=(XZ(I)+RHB)*DIP 
VR=DSQRT(VEF**2+REFF**2*WO**2 ) 

PH I =DATAN ( VEF/REFF/WO ) 

Q{ I )=- . 5D0*DC0S(PHI ) *RH0*VR**2*YI ( I ) *REFF*DIP 
800 CONTINUE 
QI=Q(ND) 

ZI=XZ(ND) 

SUMQ1=0.0D0 

CALL SPLINT(NP2,XZ,Q,SUMQ1) 

XZ(1)=0.D0 

Q(1)=QI 

DO 830 J=2,NP3 
XZ(J)=(J-1)*R2/FL0AT(NP3-1) 

JJ=ND+J-1 

YI( JJ)=2.D0*DSQRT(B3*B3*(1-(XZ( J)-B4)**2/B4**2) ) 

REFF=(XZ(J)+R1)*DIP 

VR=DSQRT( VEF**2+REFF**2*WO**2 ) 

PHI =DATAN ( VEF/REFF/WO ) 

Q( J ) =- . 25D0*DC0S ( PHI ) *RHO*VR**2*YI ( JJ ) *REFF*DIP 
830 CONTINUE 

CALL SPLINT(NP3,XZ,Q,SUMQ1) 

CTRQ=SUMQ1 
750 CONTINUE 
RETURN 
END 
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SUBROUTINE FACTRS ( SKPl , SKP2 , I S , ZETA, NP , J , XHB ,YI,YIBZ,YIAZ) 

THIS SUBROUTINE CALCULATES ALL THE INDUCTION FACTORS NEEDED 

THE INPUT PARAMETERS IN THE ARGUMENT ARE: 

ZETA RADIAL COORDINATES 

NP NUMBER OF POINTS USED IN THE SPLINE INTERPOLATION 

J INDEX REFERING TO THE CORRESPONDING POINT ON THE BLADE 

THE OUTPUT PARAMETERS OUT THE ARGUMENT ARE: 

YI INDUCTION FACTOR 

YIBZ INDUCTION FACROR AT ZETP-DEL 

YIAZ INDUCTION FACROR AT ZETP+DEL 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION SKP1(9) ,SKP2(9) ,ZETA(81) , IS(81) , AA(4, 81) ,XZ(81) 
1,Q(81),YI(81) 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 
DEL1=0.000001D0 
IF(NB.EQ.l) GO TO 10 
TT=TTP 
N0PT=2 

CALL CALCI(P2,N0PT) 

P2 IS THE INTEGRAL FROM 0 TO INFINITY OF THE INTEGRAND M 
WHERE 12= (ZETA-ZETP)* INTEGRAL OF M 
10 NOPT=l 
NAB=0 
NCTR=0 
ICTR=0 

INTERPOLATES VALUES OF I BETWEEN SKPl(J) & SKP2(J) WHERE 
SKPl(J) & SKP2(J) ARE ALREADY READ IN FOR EACH ZETP LOCATION 
DO 20 1=1, NP 
TT=ZETA(I) 

FIND THE VALUE OF DO LOOP PARAMETER FOR WHICH ZETA=ZETP 
IF(DABS(TTP-TT) .LT.O.OOOOOIDO) ICTR=I 
STORE THE INDUCTION FACTOR IN ARRAY YI 

IF(TT.LT.SKP1( J) .OR. TT.GT. SKP2 ( J) ) CALL CALCI(XI,NOPT) 

STORE WHICH INDUCTION FACTOR CALCULATIONS HAVE BEEN SKIPPED 
IF(TT.LT.SKP1(J) .OR. TT. GT . SKP2 ( J ) ) YI(I)=XI 
IF(TT.GE.SKP1( J) .AND. TT. LE . SKP2 ( J) ) IS(I)=1 
20 CONTINUE 

DO 30 1=1, NP 
ICK=0 

TT=ZETA(I) 

IF(IS(I) .NE.l .AND. TT.GE.SKP1( J)-.085D0 .AND. TT . LE . SKP2 ( J) 
1+.085D0 .OR. I.EQ.ICTR) ICK=2 

NAB IS THE NUMBER OF POINTS AT WHICH THE DIRECT CALCULATION OF 
THE INDUCTION FACTORS HAVE BEEN PERFORMED 
IF(ICK.E0.2) NAB=NAB+1 
IF(ICK.EQ.2) XZ(NAB)=TT 
C STORE THE COMPUTED INDUCTION FACTORS IN ARRAY Q 

IF(ICK.EQ.2 .AND. I.NE.ICTR) Q(NAB)=YI(I) 

IF( I .EQ. ICTR) Q(NAB)=1.0D0 
30 IF( I.EQ.ICTR) NCTR=NAB 
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USE THE SPLINE INTERPOLATION^ METHOD TO OBTAIN INDUCTION 
FACTORS THAT WERE SKIPPED IN THE DIRECT CALCULATIONS 
SPLCOE IS A COMPUTER SUBROUTINE WHICH GENERATES SPLINE 
INTERPOLATION COEFFICIENTS. 

CALL SPLCOE(XZ,Q,NAB,AA) 

DO 40 1=1, NP 
TT=ZETA( I) 

USE SPLINE INTERPOLATION FORMULA 
IF(IS(I) .EQ.l .AND. I.LT.ICTR) YI(I)= 

1 AA(1,NCTR-1)*(XZ(NCTR)-TT)**3 
2+AA(2,NCTR-l)*(TT-XZ(NCTR-l) )**3 
3+AA(3,NCTR-l)*(XZ(NCTR)-TT) 

4+AA ( 4 , NCTR- 1 ) * ( TT-XZ ( NCTR- 1 ) ) 

IF(I.EQ.ICTR) YI(I)=1.0D0 
40 IF(IS(I) .EQ.l .AND. I.GT.ICTR) 

1 YI(I)=AA(1,NCTR)*(XZ(NCTR+1)-TT)**3 

2 +AA(2,NCTR)*(TT-XZ(NCTR) )**3 

3 +AA(3,NCTR)*(XZ(NCTR+1)-TT) 

4 +AA( 4, NCTR)* (TT-XZ (NCTR) ) 

IF(TTP.GT.XHB+DEL1) TT=TTP-DEL1 
I F ( TTP . GT . XHB+DELl ) 

1 Y I BZ=AA ( 1 , NCTR- 1 ) * ( XZ ( NCTR ) - TT ) * * 3 

2 +AA(2,NCTR-1)*(TT-XZ(NCTR-1) )**3 

3 +AA(3,NCTR-1)*(XZ(NCTR)-TT) 

4 +AA(4,NCTR-1)*(TT-XZ(NCTR-1) ) 

I F ( 1 . ODO-TTP . GT . DELI ) TT=TTP+DEL1 

I F ( 1 . ODO-TTP . GT . DELI ) 

1 Y I AZ=AA ( 1 , NCTR ) * ( XZ ( NCTR+ 1 ) -TT ) * * 3 

2 +AA(2,NCTR)*(TT-XZ(NCTR) )**3 

3 +AA(3,NCTR)*(XZ(NCTR+1)-TT) 

4 +AA( 4, NCTR)* (TT-XZ (NCTR) ) 

ALL INDUCTION FACTORS HAVE NOW BEEN CALCULATED AND STORED 
IN YI, YIBZ, YIAZ. YIBZ & YIAZ ARE THE INDUCTION FACTORS 
AT ZETP-DEL AND ZETP+DEL RESPECTIVELY. 

RETURN 

END 
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SUBROUTINE SPLINT ( ND, ZZ, C, SUM) 

THIS SUBROUTINE COMPUTES A FINITE INTEGRAL 
IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION ZZ(81) ,ZC(4,81) ,C(81) 

CALL SPLCOE(ZZ,C,ND,ZC) 

NN=ND-1 
SUM=O.DO 
DO 10 K=1,NN 

10 SUM=SUM+.25D0*ZC(1,K)*(ZZ(K+1)-ZZ(K) )**4 

1 +.25D0*ZC(2,K)*(ZZ(K+1)-ZZ(K) )**4 

2 +.50D0*ZC(3,K)*(ZZ(K+1)-ZZ(K) )**2 

3 +.50D0*ZC(4,K)*(ZZ(K+1)-ZZ(K) )**2 
WRITE (7, 100) SUM,NN 

WRITE(7,200) (ZZ(I),C(I),I=1,NN) 

100 FORMAT(F10.4, 110) 

200 FORMAT (4F10. 4) 

RETURN 

END 
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SUBROUTINE TINGRL ( NS, NP , ZETA,XHB, DELI, YI , YIBZ , YIAZ, TINT) 
IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/CCC/P2 , TT , TTP , P I , EL , NB 

DIMENSION ZETA(81),S(81),XZ(81),Q(81),YI(81) 

DEL=DEL1 

ND=0 

ARG=NS*P I * ( TTP-XHB ) / ( 1 . DO-XHB ) 

IF(TTP.LE.DEL+XHB) GO TO 30 

THIS SECTION CALCULATES THE INTEGRAL FROM XHB TO ZETP-DEL 

DO 10 1=1, NP 
ND=ND+1 

IF(ZETA(I)-TTP.GE.-DEL1) GO TO 20 
XZ(I)=ZETA(I) 

ARGP= ( ZETA( I ) -XHB ) / ( 1 . ODO-XHB ) 

STORE THE INTEGRAND IN ARRAY Q 
10 Q ( I ) =DCOS (NS*P I *ARGP ) *YI ( I ) /( ZETA ( I ) -TTP ) 

20 XZ(ND)=TTP-DEL 

ARGP= ( XZ ( ND ) -XHB ) / ( 1 . ODO -XHB ) 
Q(ND)=DCOS(NS*PI*ARGP)*YIBZ/(XZ(ND)-TTP) 

INTEGRATE Q USING SPLINE INTEGRATION FORMULA 
SUMQ=0 . ODO 

CALL SPLINT (ND,XZ,Q,SUMQ) 

30 CONTINUE 

IF( TTP-XHB. LT. DEL .OR. 1 . ODO-TTP . LT.DEL) DEL=DEL/2 . ODO 


THIS SECTION CALCULATES THE INTEGRAL FROM ZETP-DEL TO ZETP+DEL 

SUMR=0 . ODO 
ATERM=O.DO 
BTERM=O.DO 
IF(NB.EQ.l) GO TO 40 
ATERM=DCOS( ARG) *2 . 0D0*DEL*P2 
40 BTERM=-NS*PI*2 . ODO*DEL*DSIN( ARG) 

SUMR=ATERM+BTERM 
DEL=DEL1 

IF(1. ODO-TTP. LE. DEL) GO TO 60 

THIS SECTION CALCULATES THE INTEGRAL FROM ZETP+DEL TO 1 
ND=1 

XZ(ND)=TTP+DEL 

ARGP= ( XZ ( ND ) -XHB )/( 1 . ODO-XHB ) 
S(ND)=DCOS(NS*PI*ARGP)*YIAZ/(XZ(ND)-TTP) 

DO 50 1=1, NP 

IF (ZETA(I)-TTP.LE.DELl) GO TO 50 
ND=ND+1 
XZ(ND)=ZETA(I) 

C STORE THE INTEGRAND IN ARRAY S 
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ARGP= ( XZ ( ND ) -XHB ) / ( 1 . ODO-XHB ) 
S(ND)=DCOS(NS*PI*ARGP)*YI(I)/(XZ(ND)-TTP) 

50 CONTINUE 

C INTEGRATE S USING SPLINE INTEGRATION FORMULA 
SUMS=O.ODO 

CALL SPLINT(ND,XZ,S,SUMS) 

C TINT IS THE TOTAL INTEGRAL FROM XHB TO 1 
60 TINT=SUMQ+SUMR+SUMS 
RETURN 
END 
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n o o 


SUBROUT INK SK I P ( ZETP , RB , OMEGA ,VEL,N,SKP1,SKP2) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION ZETP(9) , SKP1(9) , SKP2(9) 

WRITE(7,501) VEL 

501 FORMAT( ' HERHE VEL ',F10.4) 

WRITE(7,502) ( ZETP ( I ) , 1=1 , N) 

502 FORMAT( ' ZETAP VAL ',F10.4) 

SKP1(1)=ZETP(1) 

SKP2(N)=1.0D0 

FR=RB*OMEGA/VEL*0 . O020973D0 
IF(FR.GT.l.) FR=1.0 
DO 50 J=2,N 

SKP1( J)=ZETP(J) -0.050-0. 025*ZETP(J)**2-.025*ZETP(J)**2*FR**3 
IF(SKP1( J) .LT.ZETP(l) ) SKP1(J)= ZETP ( 1 ) +0 . 0001 
50 CONTINUE 

DO 100 J=1,N-1 

SKP2( J)=ZETP( J)+0.050+.025*(ZETP(N-l) )**2 
1 + . 025 * ( ZETP ( J ) /ZETP ( N- 1 ) ) * *2 * FR* * 3 

IF(SKP2( J) .GT.ZETP(N) ) SKP2(J)= ZETP (N) -0 . 0001 
100 CONTINUE 

WRITE (7,200) ( SKPl ( I ) , SKP2 ( I ) , 1=1 , N ) 

200 FORMAT (2F12. 6) 

RETURN 

END 
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